本帖最后由 大禹man 于 2017-2-25 16:07 编辑

需要拟合的函数为zeta，用omeg的前两个数计算得到系数a0和a1，他们会在循环过程中随着miu0和miu1变化，我的程序如下，可以运行但是无法得到真实的迭代值，最终的迭代值应该为miu1=500,miu2=100,请大神帮忙看看:handshake

%算例

%*************************************************************

%状态空间模态叠加法

%=============================================================

clc

clear

clear all

freedom=6;%自由度数

m=1;

k=40000;

alpha=20;beta=0.00003;

miu1=500;miu2=100;              %松弛因子

I=eye(freedom);

zero=zeros(freedom);

M=m*I;

K=k*(diag(repmat([2], 1, freedom))+diag(repmat([-1], 1, freedom-1), 1)+diag(repmat([-1], 1, freedom-1), -1));

C1=alpha*M;C2=beta*K;        %阻尼系数矩阵

% ======写为状态空间形式=====

P=[C1+C2 M -C1/miu1 -C2/miu2;M zero zero zero;-C1/miu1 zero C1/miu1^2 zero;-C2/miu2 zero zero C2/miu2^2];

Q=-[-K zero zero zero;zero M zero zero;zero zero -C1/miu1 zero;zero zero zero -C2/miu2];

%=====求解复模态====

LQ=length(Q);

[Y,Lambda]=eig(-inv(P)*Q);

Lambda=diag(Lambda);

temp=Lambda(find(imag(Lambda)>0));

omeg1=abs(temp)/2/pi;%得到无阻尼固有频率

zeta1=-real(temp)./abs(temp)*100;%阻尼比

omeg=sort(omeg1);omeg=omeg.';

zeta=sort(zeta1);zeta=zeta.';

save zetaomeg zeta omeg

load zetaomeg

oomeg1=omeg(:,1:2);%前两列用来计算比例系数

oomeg2=omeg(:,3:5);%后四列用来迭代

zzeta1=zeta(:,1:2);%前两列用来计算比例系数

zzeta2=zeta(:,3:5);%后四列用来迭代

% 计算函数f的雅克比矩阵，是解析式

syms miu11 miu22

temp=zzeta1.';%前两阶阻尼比

temp1=1/(1+(oomeg1(:,1)/miu11)^2);temp2=1/(1+(oomeg1(:,1)/miu22)^2);

temp3=1/(1+(oomeg1(:,2)/miu11)^2);temp4=1/(1+(oomeg1(:,2)/miu22)^2);

TEMP=[temp1/oomeg1(:,1) omeg(:,1)*temp2;temp3/oomeg1(:,2) omeg(:,2)*temp4];

TTEMP=2*inv(TEMP)*temp;Alpha=TTEMP(1);Beta=TTEMP(2);    %比例系数

%写出误差项(总共有四项)

r1=Alpha/(1+(oomeg2(:,1)/miu11)^2)/2/oomeg2(:,1)+Beta*oomeg2(:,1)/(1+(oomeg2(:,1)/miu22)^2)/2-zzeta2(:,1);

r2=Alpha/(1+(oomeg2(:,2)/miu11)^2)/2/oomeg2(:,2)+Beta*oomeg2(:,2)/(1+(oomeg2(:,2)/miu22)^2)/2-zzeta2(:,2);

r3=Alpha/(1+(oomeg2(:,3)/miu11)^2)/2/oomeg2(:,3)+Beta*oomeg2(:,3)/(1+(oomeg2(:,3)/miu22)^2)/2-zzeta2(:,3);

%写出阻尼比解析式

z1=r1+zzeta2(:,1);z2=r2+zzeta2(:,2);z3=r3+zzeta2(:,3);

J1=jacobian(z1,[miu11 miu22]);

J2=jacobian(z2,[miu11 miu22]);

J3=jacobian(z3,[miu11 miu22]);

% 2. LM算法

% 初始猜测s

miu10=10000; miu20=5000;

r1_init=double(subs(r1,[miu11,miu22],[miu10,miu20]));

r2_init=double(subs(r2,[miu11,miu22],[miu10,miu20]));

r3_init=double(subs(r3,[miu11,miu22],[miu10,miu20]));

y_init=[r1_init r2_init r3_init];

% 数据个数

Ndata=length(zzeta2);

% 参数维数

Nparams=2;

% 迭代最大次数

n_iters=50;

% LM算法的阻尼系数初值

lamda=0.001;

% step1: 变量赋值

updateJ=1;

miu11_est=miu10;

miu22_est=miu20;

% step2: 迭代

for it=1:n_iters

if updateJ==1

% 根据当前估计值，计算雅克比矩阵

J11=double(subs(J1,[miu11,miu22],[miu11_est,miu22_est]));

J22=double(subs(J2,[miu11,miu22],[miu11_est,miu22_est]));

J33=double(subs(J3,[miu11,miu22],[miu11_est,miu22_est]));

J=[J11;J22;J33];

% 根据当前参数，得到函数值

r1_est=double(subs(r1,[miu11,miu22],[miu11_est,miu22_est]));

r2_est=double(subs(r2,[miu11,miu22],[miu11_est,miu22_est]));

r3_est=double(subs(r3,[miu11,miu22],[miu11_est,miu22_est]));

%得到误差

y_est=[r1_est r2_est r3_est];

% 计算(拟)海塞矩阵

H=J'*J;

% 若是第一次迭代，计算误差

if it==1

e=dot(y_est,y_est);

end

end

% 根据阻尼系数lamda混合得到H矩阵

H_lm=H+(lamda*eye(Nparams,Nparams));

% 计算步长dp，并根据步长计算新的可能的\参数估计值

g=-J.'*y_est(:);

dp=(H_lm)\g;

miu11_lm=miu11_est+dp(1);

miu22_lm=miu22_est+dp(2);

if norm(dp)<1e-5

break

end

% 计算新的可能估计值对应的zeta和计算残差e

r1_est_lm=double(subs(r1,[miu11,miu22],[miu11_lm,miu22_lm]));

r2_est_lm=double(subs(r2,[miu11,miu22],[miu11_lm,miu22_lm]));

r3_est_lm=double(subs(r3,[miu11,miu22],[miu11_lm,miu22_lm]));

y_est_lm=[r1_est_lm r2_est_lm r3_est_lm];

e_lm=dot(y_est_lm,y_est_lm);

% 根据误差，决定如何更新参数和阻尼系数

%     while abs(miu11_lm-miu11_est)<10|abs(miu22_lm-miu22_est)<10

%         break

%     end

if e_lm

lamda=lamda/10;

miu11_est=miu11_lm;

miu22_est=miu22_lm;

e=e_lm;

disp(e);

updateJ=1;

else

updateJ=0;

lamda=lamda*10;

end

end

%显示优化的结果

miu11_est

miu22_est

2017-2-25 16:07 上传

其中miu0和miu1为需要拟合的元素，Omega和zeta为自变量和应变量，共有三组数据 ... ...

2.png

(24.34 KB, 下载次数: 3)

2017-2-25 16:07 上传

系数a0和a1的计算公式