Matrix close to singular or badly scaled problem in Gauss-Newton LMS calculations
    3 次查看(过去 30 天)
  
       显示 更早的评论
    
I am trying to implement a function that minimizes the square of the residuals (LMS) for a nonlinear function. I am trying to use the Gauss-Newton method. The nonlinear function is the Fresnel reflection coeficients, which are calculated with another function I implemented, called bwTMM.
The minimizing function is called GNTMM. When I try the function with several different parameters i get the error: 
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
> In GNTMMs (line 58)
Any ideas? I insert 3 codes: The function bwTMM that calculates the Fresnel coefficients, the function that tries to fit the parameters, and the function in which I try the fitting function and get an error.
Thank you for taking the time to read:
function [GMs,GMp,phiout,Rs,Ts,Rp,Tp]=bwTMM(n,e,phi0,lambda)%e y lambda en mismas unidades
ko=2*pi/lambda;%calculo ko
phi=deg2rad(phi0);%defino el primer phi (angulo de incidencia)
Ms=eye(2);
Mp=eye(2);
for k=1:length(e)
    phi=asin((n(k)/n(k+1))*sin(phi));
    p=n(k+1)*cos(phi);%p es sqrt(epsilon/mu)cos(phi), pero si medio no magnetico mu=1 y, como sqrt(epsilon)=n, queda p
    B=ko*n(k+1)*e(1)*cos(phi);
    q=(1/n(k+1))*cos(phi);% Esto es para TM, o polarizacion paralela al plano de incidencia
    %--------Matriz para polarizacion TE (plarizacion s)
    ms11=cos(B);
    ms12=-(1i/p)*sin(B);
    ms21=-1i*p*sin(B);
    ms22=cos(B);
    Ms=Ms*[ms11,ms12;ms21,ms22];
    %------Matriz polarizacion TM (p polarized)-------
    mp11=cos(B);
    mp12=-(1i/q)*sin(B);
    mp21=-1i*q*sin(B);
    mp22=cos(B);
    Mp=Mp*[mp11,mp12;mp21,mp22];
end
phiout=rad2deg(asin((n(length(e)+1)/n(length(n)))*sin(phi)));
p1=n(1)*cos(phi0);
pl=n(end)*cos(phiout);
q1=(1/n(1))*cos(phi0);
ql=(1/n(end))*cos(phiout);
%-----Extraigo elementos de Ms global:-------
ms11g=Ms(1,1);
ms12g=Ms(1,2);
ms21g=Ms(2,1);
ms22g=Ms(2,2);
%------Extraigo elementos de Mp global:------
mp11g=Mp(1,1);
mp12g=Mp(1,2);
mp21g=Mp(2,1);
mp22g=Mp(2,2);
%----Coeficientes R Y T para polarizacion s------
%-----Calculo de r y t (coeficientes de amplitud)
dens=(ms11g+pl*ms12g)*p1+(ms21g+pl*ms22g);
rs=((ms11g+pl*ms12g)*p1-(ms21g+pl*ms22g))/dens;
ts=2*p1/dens;
%-----Coeficientes R y T (intensidad)-------
Rs=(abs(rs))^2;
Ts=(pl/p1)*(abs(ts))^2;
%----Coeficientes R Y T para polarizacion p------
%-----Calculo de r y t (coeficientes de amplitud)
denp=(mp11g+ql*mp12g)*q1+(mp21g+ql*mp22g);
rp=((mp11g+ql*mp12g)*q1-(mp21g+ql*mp22g))/denp;
tp=2*q1/denp;
%-----Coeficientes R y T (intensidad)-------
Rp=(abs(rp))^2;
Tp=(ql/q1)*(abs(tp))^2;
%------ Matrices ------
GMs=Ms;
GMp=Mp;
%----- Code 2, fitting function------------
function [parameters]=GNTMMs(data_points,nvalues,evalues,lambda,parameters,h,max_counts,min_step)
counts=0;
found=false;
angle=data_points(:,1);
Rsvalues=data_points(:,2);
while (found==false & counts<max_counts)
    %------Calculo de los coeficientes para estos parametros------
    for k=1:length(angle)
    phi0=angle(k);
    [~,~,~,Rs,~,Rp,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3),evalues(1)],phi0,lambda);
    Rfp(k)=Rp;
    Rfs(k)=Rs;
    %Elementos del jacobiano---
    [~,~,~,Rs1p,~,Rp1p,~]=bwTMM([1,parameters(1)+h(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3),evalues(1)],phi0,lambda);
    [~,~,~,Rs1m,~,Rp1m,~]=bwTMM([1,parameters(1)-h(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3),evalues(1)],phi0,lambda);
    Rfp1p(k)=Rp1p;
    Rfs1p(k)=Rs1p;
    Rfp1m(k)=Rp1m;
    Rfs1m(k)=Rs1m;
    %---------------------------------
    [~,~,~,Rs2p,~,Rp2p,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2)+h(2),parameters(3),evalues(1)],phi0,lambda);
    [~,~,~,Rs2m,~,Rp2m,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2)-h(2),parameters(3),evalues(1)],phi0,lambda);
    Rfp2p(k)=Rp2p;
    Rfs2p(k)=Rs2p;
    Rfp2m(k)=Rp2m;
    Rfs2m(k)=Rs2m;
    %------------------------
    [~,~,~,Rs3p,~,Rp3p,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3)+h(2),evalues(1)],phi0,lambda);
    [~,~,~,Rs3m,~,Rp3m,~]=bwTMM([1,parameters(1),nvalues(1),nvalues(2),1],[parameters(2),parameters(3)-h(2),evalues(1)],phi0,lambda);
    Rfp3p(k)=Rp3p;
    Rfs3p(k)=Rs3p;
    Rfp3m(k)=Rp3m;
    Rfs3m(k)=Rs3m;
    end
    %------Calculo jacobiano
    for r=1:length(Rsvalues)
        J(r,1)=(Rfs1p(r)-Rfs1m(r))/(2*h(1));
        J(r,2)=(Rfs2p(r)-Rfs2m(r))/(2*h(2));
        J(r,3)=(Rfs3p(r)-Rfs3m(r))/(2*h(2));
    end
     residuals=(Rsvalues-transpose(Rfs))*1000000;
     Jt=transpose(J);
     stepgn=(1000000*(Jt*J))\(-Jt*residuals);
     parameters=parameters+transpose(stepgn)/1000000;
     if norm(stepgn)<min_step
         found=true;
     end
     counts=counts+1;
end
end
%------ Code 3--------
close all
clear
clc
data=load('4mM0cs1data.txt');
angle=data(:,1);
Exrp=data(:,2);
Exrs=data(:,3);
Exratio=data(:,4);
%-----known values-----
nvidrio=1.5116+1i*(9.1316)*10^(-9);
evidrio=(1.25)*10^(-6);%nm
nFTO=1.8;
%----------
[parameters]=GNTMMs([angle,Exrs],[nFTO,nvidrio],[evidrio],760,[3,500,200],[0.4,200],50,0.01);
1 个评论
  Torsten
      
      
 2022-3-16
				Why don't you use one of the MATLAB routines for the least squares part, e.g. lsqnonlin ?
回答(0 个)
另请参阅
类别
				在 Help Center 和 File Exchange 中查找有关 Guidance, Navigation, and Control (GNC) 的更多信息
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!

