24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 584  |  回复: 1

我爱信息光学

新虫 (初入文坛)

[求助] 最小二乘法编程遇到问题求帮助

我用阻尼最小二乘法求解非线性方程组f=acos(d-b)+c最优值问题,d是自变量,结果我想得到迭代次数和a,b,c在迭代中的变化,该怎么修改,求大神帮忙,代码如下:
clear all;
tic
vx=zeros(1,90);
vy=vx;
vz=vx;
%%初始值
vx(1,=8;
vy(1,=8;
vz(1,=1;
%%加噪声
vxn=awgn(vx(1,,0);
vyn=awgn(vy(1,,0);
vzn=awgn(vz(1,,0);
%vxn=3*6^(0.5)+randn(90,1);
%vyn=pi/4+randn(90,1);
%vzn=0.5+randn(90,1);
a=((vxn).^2.+(vyn).^2).^(0.5)*3^(0.5)*0.5;
b=atan(vyn./vxn);
c=0.5*vzn;
%%方位角,90个点
dd=pi/45:pi/45:2*pi;
vr=zeros(1,90);
vro=zeros(1,90);
for i=1:1:90
    vr(i)=a(i)*cos(dd(i)-b(i))+c(i);
    vro(i)=4*6^(0.5)*cos(dd(i)-pi/4)+0.5;
end
%%
% 计算函数f的雅克比矩阵,是解析式
syms a b c y x real;
y=a*cos(x-b)+c;
Jsym=jacobian(y,[a b c])
% 拟合用数据。参见《数学试验》,p190,例2
data_1=dd;
obs_1=vr;
% 2. LM算法
% 初始
%a0=3*6^(0.5);
%b0=pi/4;
%c0=0.5;
a0=1;
b0=8;
c0=0;
y_init = a0*cos(data_1(1)-b0)+c0;
% 数据个数
Ndata=length(obs_1);
% 参数维数
Nparams=3;
% 迭代最大次数
n_iters=50;
% LM算法的阻尼系数初值
lamda=0.1;
% step1: 变量赋值
updateJ=1;
a_est=a0;
b_est=b0;
c_est=c0;
% step2: 迭代
for it=1:n_iters
    if updateJ==1
        % 根据当前估计值,计算雅克比矩阵
        J=zeros(Ndata,Nparams);
        for i=1:length(data_1)
            %J(i,=[exp(-b_est*data_1(i)) -a_est*data_1(i)*exp(-b_est*data_1(i))];
            J(i,=[cos(data_1(i)-b_est) a_est*sin(data_1(i)-b_est) 1];
        end
        % 根据当前参数,得到函数值
        %y_est = a_est*exp(-b_est*data_1);
        y_est=a_est*cos(data_1-b_est)+c_est;
        % 计算误差
        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(;
    a_lm=a_est+dp(1);
    b_lm=b_est+dp(2);
    c_lm=c_est+dp(3);
    % 计算新的可能估计值对应的y和计算残差e
    %y_est_lm = a_lm*exp(-b_lm*data_1);
    y_est_lm=a_lm*cos(data_1-b_lm)+c_lm;
    d_lm=obs_1-y_est_lm;

    e_lm=dot(d_lm,d_lm);
    % 根据误差,决定如何更新参数和阻尼系数
    if e_lm<e
        lamda=lamda/10;
        a_est=a_lm;
        b_est=b_lm;
        c_est=c_lm;
        e=e_lm;
        it=it+1
        disp(e);
        updateJ=1;
    else
        updateJ=0;
        lamda=lamda*10;

    end
end
dis=[a_est/((1+tan(b_est)^2)^(0.5)*3^(0.5)/2) a_est*tan(b_est)/((1+tan(b_est)^2)^(0.5)*3^(0.5)/2) c_est*2]
toc

%%

vrn=a_est*cos(dd-b_est)+c_est;
plot(dd,vr,'r+');
hold on;
plot(dd,vrn);
plot(dd,vro,'g');
hold off;
回复此楼

» 猜你喜欢

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

manmanbobo

新虫 (初入文坛)

应该不难的,你就改个代码就行了。
2楼2017-10-11 09:19:21
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 我爱信息光学 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 070305高分子化学与物理 304分求调剂 +7 c297914 2026-03-28 7/350 2026-03-29 01:17 by 544594351
[考研] 081200-11408-276学硕求调剂 +6 崔wj 2026-03-26 6/300 2026-03-29 01:11 by hanserlol
[考研] 070300化学354求调剂 +3 101次希望 2026-03-28 3/150 2026-03-29 00:28 by 544594351
[考研] 11408软件工程求调剂 +3 Qiu学ing 2026-03-28 3/150 2026-03-28 21:50 by zhq0425
[考研] 332求调剂 +4 @MZB382400 2026-03-28 4/200 2026-03-28 21:02 by 唐沐儿
[考研] 343求调剂 +5 爱羁绊 2026-03-28 5/250 2026-03-28 20:53 by 唐沐儿
[考研] 347求调剂 +3 山顶见α 2026-03-25 3/150 2026-03-28 14:13 by 唐沐儿
[考研] 291求调剂 +6 HanBeiNingZC 2026-03-24 6/300 2026-03-28 07:55 by baoball
[考研] 307求调剂 +8 超级伊昂大王 2026-03-24 9/450 2026-03-27 15:34 by 超级伊昂大王
[考研] 复试调剂,一志愿南农083200食品科学与工程 +5 XQTJZ 2026-03-26 5/250 2026-03-27 14:49 by 狂炫麦当当
[考研] 085600,材料与化工321分调剂 +4 大馋小子 2026-03-27 6/300 2026-03-27 14:11 by 松花缸1201
[考研] 0703化学338求调剂! +6 Zuhui0306 2026-03-26 7/350 2026-03-27 10:35 by shangxh
[考研] 294分080500材料科学与工程求调剂 +4 柳溪边 2026-03-26 4/200 2026-03-26 21:14 by XPU李庆
[考研] 321求调剂 +6 Ymlll 2026-03-24 6/300 2026-03-26 20:50 by 不吃魚的貓
[考研] 085601求调剂总分293英一数二 +4 钢铁大炮 2026-03-24 4/200 2026-03-26 16:28 by dick_runner
[考研] 297求调剂 +6 田洪有 2026-03-26 6/300 2026-03-26 15:55 by 不吃魚的貓
[考研] 材料与化工304求B区调剂 +3 邱gl 2026-03-25 3/150 2026-03-25 19:03 by Ainin_
[考研] 生物技术与工程 +3 1294608413 2026-03-25 4/200 2026-03-25 18:02 by 1294608413
[考研] B区考研调剂 +4 yqdszhdap- 2026-03-22 5/250 2026-03-25 08:51 by baoball
[考研] 300分,材料,求调剂,英一数二 +5 超赞的 2026-03-24 5/250 2026-03-24 21:07 by 星空星月
信息提示
请填处理意见