In solving two non-linear equations using fsolve, I am getting an error 'Failure in Initial objective function evaluation. FSOLVE cannot continue'. How to resolve?

%The code below solves two equations. The initial part of the code generates 10000 values of random data U,V,x_0,y_0 which are seed values for fsolve


clc
clear all
close all
n=10000;
%Generating 10000 values of Gaussian distributed random delay in the range
%5*10^-3 and 6*10^-3
a=5*10^-3;
b=6*10^-3;
X = (b-a).*randn(1,n) + a;

X1 = X;
d=10*10^-3;

%Generating T_1_i_A(:,1)
a=0.005;
b=5;
x = (b-a).*rand(1,n)+a;

T_1_i_A=zeros(n,30);
T_1_i_A(:,1)=transpose(x);

%Generating T_1_i_A for N=2:30 where N is no. of rounds of message exchange
for j=2:30
    T_1_i_A(:,j)=T_1_i_A(:,j-1)+X(j);
end

%Generating T_2_i_B and T_2_i_P from T_1_i_A,d,X
T_2_i_B=zeros(n,30);

for j=1:30
 T_2_i_B(:,j)=T_1_i_A(:,j)+20*10^-6+1.01*(T_1_i_A(:,j)-T_1_i_A(:,1))+15*10^-3+transpose(X);
end

T_2_i_P=zeros(n,30);

for j=1:30
 T_2_i_P(:,j)=T_1_i_A(:,j)+30*10^-6+1.105*(T_1_i_A(:,j)-T_1_i_A(:,1))+5*10^-3+transpose(X1);
end

%finding U and V
U=T_2_i_P-T_2_i_B;
V=T_1_i_A-T_1_i_A(:,1);

fun = @root2d;

%Generating 10000 random values of pre-defined clock offset in the range
%5*10^-6 and 30*10^-6
a=5*10^-6;
b=30*10^-6;
n=10000;
x_0 = (b-a).*rand(1,n) + a;

%Generating 10000 random values of pre-defined clock skew in the range
%1.01 and 1.20
a1=1.01;
b1=1.20;
y_0=(b1-a1).*rand(1,n)+a1;


x0=zeros(2,n);
x=zeros(2,n);
for i=1:n
    x0(:,i)=[x_0(i); y_0(i)];
end
for N=1:30
n=10000;
for i=1:n
    x(:,i)= fsolve(@root2d,x0(:,i),[],U(i,N),V(i,N),N);
end
end



function f = root2d(x)

f(1)=(1.5./(U-x(1)-x(2)*V+10^-4)-0.05./(U-x(1)-x(2)*V+10^-4).^2)*N+0.1*sum(1:N);

f(2)=((1.5*V)./(U-x(1)-x(2)*+10^-4)-0.05*V./(U-x(1)-x(2)*V+10^-4).^2+0.1*V)*N;
return 
end

Answer 1

[Walter Roberson]

You really need to recheck your root2d code. You are passing scalar U and V entries but your for i loop involves U(:,j) which requires non-scalar U. But if U is non-scalar then your accumulated k is going to be non-scalar and then storing it into f(1) and f(2) would fail.

Your for j loop accumulating into k uses expressions that are independent of j except for the final 0.1*j -- so you might as well only calculate the expression once and then add on the cumulative total sum(1:N) * 0.1 .

Your for i loop is independent of i so it is going to give a result that is N times the value of the expression once.

Note that your use of U(:,j) inside for i is not indexing by i and is instead using whatever value of j happens to be sitting around in the workspace -- in particular it would be using U(:,N) since for j=1:N is going to leave j set to N after the loop.

%The code below solves two equations. The initial part of the code generates 10000 values of random data U,V,x_0,y_0 which are seed values for fsolve



n=10000;
%Generating 10000 values of Gaussian distributed random delay in the range
%5*10^-3 and 6*10^-3
a=5*10^-3;
b=6*10^-3;
X = (b-a).*randn(1,n) + a;

X1 = X;
d=10*10^-3;

%Generating T_1_i_A(:,1)
a=0.005;
b=5;
x = (b-a).*rand(1,n)+a;

T_1_i_A=zeros(n,30);
T_1_i_A(:,1)=transpose(x);

%Generating T_1_i_A for N=2:30 where N is no. of rounds of message exchange
for j=2:30
    T_1_i_A(:,j)=T_1_i_A(:,j-1)+X(j);
end

%Generating T_2_i_B and T_2_i_P from T_1_i_A,d,X
T_2_i_B=zeros(n,30);

for j=1:30
 T_2_i_B(:,j)=T_1_i_A(:,j)+20*10^-6+1.01*(T_1_i_A(:,j)-T_1_i_A(:,1))+15*10^-3+transpose(X);
end

T_2_i_P=zeros(n,30);

for j=1:30
 T_2_i_P(:,j)=T_1_i_A(:,j)+30*10^-6+1.105*(T_1_i_A(:,j)-T_1_i_A(:,1))+5*10^-3+transpose(X1);
end

%finding U and V
U=T_2_i_P-T_2_i_B;
V=T_1_i_A-T_1_i_A(:,1);

fun = @root2d;

%Generating 10000 random values of pre-defined clock offset in the range
%5*10^-6 and 30*10^-6
a=5*10^-6;
b=30*10^-6;
n=10000;
x_0 = (b-a).*rand(1,n) + a;

%Generating 10000 random values of pre-defined clock skew in the range
%1.01 and 1.20
a1=1.01;
b1=1.20;
y_0=(b1-a1).*rand(1,n)+a1;


x0=zeros(2,n);
x=zeros(2,n);
for i=1:n
    x0(:,i)=[x_0(i); y_0(i)];
end

options = optimoptions(@fsolve, 'display', 'none');
for N=1:30
n=10000;
for i=1:n
    x(:,i,N)= fsolve(@(x)root2d(x, U(i,N),V(i,N),N), x0(:,i), options);
end
end

function f = root2d(x, U, V, N)
k=0;
for j=1:N
    k=k+1.5./(U-x(1)-x(2)*V+10^-4)-0.05./(U-x(1)-x(2)*V+10^-4).^2+0.1*j;
end
f(1)=k;
k=0;
for i=1:N
    k=k+(1.5*V)./(U-x(1)-x(2)*+10^-4)-0.05*V./(U-x(1)-x(2)*V+10^-4).^2+0.1*V;
end
f(2)=k;
return 
end

[Siva Rama Krishna Ayaluru]

Thank you for your answer. Sorry for usage of non scalar j in i loop. I have modified the code of root2d according to your notices as below. But getting the same error.

function f = root2d(x)

f(1)=(1.5./(U-x(1)-x(2)*V+10^-4)-0.05./(U-x(1)-x(2)*V+10^-4).^2)*N+0.1*sum(1:N);

f(2)=((1.5*V)./(U-x(1)-x(2)*+10^-4)-0.05*V./(U-x(1)-x(2)*V+10^-4).^2+0.1*V)*N;
return 
end

[Walter Roberson]

I am running the code I posted. I have not received any error message yet, but it is taking a while.

Notice by the way the assignment of the results to x(:,i,N) -- so x is going to end up a 2 by 10000 x 30 array. It is not obvious to me what you are going to do with the results.

[Walter Roberson]

Minor change to improve performance by pre-allocating the output array.

I ran the previous version without error.

%The code below solves two equations. The initial part of the code generates 10000 values of random data U,V,x_0,y_0 which are seed values for fsolve



n=10000;
%Generating 10000 values of Gaussian distributed random delay in the range
%5*10^-3 and 6*10^-3
a=5*10^-3;
b=6*10^-3;
X = (b-a).*randn(1,n) + a;

X1 = X;
d=10*10^-3;

%Generating T_1_i_A(:,1)
a=0.005;
b=5;
x = (b-a).*rand(1,n)+a;

T_1_i_A=zeros(n,30);
T_1_i_A(:,1)=transpose(x);

%Generating T_1_i_A for N=2:30 where N is no. of rounds of message exchange
for j=2:30
    T_1_i_A(:,j)=T_1_i_A(:,j-1)+X(j);
end

%Generating T_2_i_B and T_2_i_P from T_1_i_A,d,X
T_2_i_B=zeros(n,30);

for j=1:30
 T_2_i_B(:,j)=T_1_i_A(:,j)+20*10^-6+1.01*(T_1_i_A(:,j)-T_1_i_A(:,1))+15*10^-3+transpose(X);
end

T_2_i_P=zeros(n,30);

for j=1:30
 T_2_i_P(:,j)=T_1_i_A(:,j)+30*10^-6+1.105*(T_1_i_A(:,j)-T_1_i_A(:,1))+5*10^-3+transpose(X1);
end

%finding U and V
U=T_2_i_P-T_2_i_B;
V=T_1_i_A-T_1_i_A(:,1);

fun = @root2d;

%Generating 10000 random values of pre-defined clock offset in the range
%5*10^-6 and 30*10^-6
a=5*10^-6;
b=30*10^-6;
n=10000;
x_0 = (b-a).*rand(1,n) + a;

%Generating 10000 random values of pre-defined clock skew in the range
%1.01 and 1.20
a1=1.01;
b1=1.20;
y_0=(b1-a1).*rand(1,n)+a1;


x0=zeros(2,n);
x=zeros(2,n);
for i=1:n
    x0(:,i)=[x_0(i); y_0(i)];
end

options = optimoptions(@fsolve, 'display', 'none');
maxN = 30;
n=10000;

x = zeros(size(x0,1), n, maxN);

for N=1:maxN
    for i=1:n
        x(:,i,N)= fsolve(@(x)root2d(x, U(i,N),V(i,N),N), x0(:,i), options);
    end
end

function f = root2d(x, U, V, N)
k=0;
for j=1:N
    k=k+1.5./(U-x(1)-x(2)*V+10^-4)-0.05./(U-x(1)-x(2)*V+10^-4).^2+0.1*j;
end
f(1)=k;
k=0;
for i=1:N
    k=k+(1.5*V)./(U-x(1)-x(2)*+10^-4)-0.05*V./(U-x(1)-x(2)*V+10^-4).^2+0.1*V;
end
f(2)=k;
return 
end

[Siva Rama Krishna Ayaluru]

Thank you very much. I want to find the Mean square error of estimated clock offset and clock skew values across the pre-defined clock offset and skew values as shown in the code below. But the MSE is fluctuating, I mean not decreasing steadily. So I tried the same program with replacing fsolve by lsqnonlin at two places. But I am getting the error 'LSQNONLIN requires the following inputs to be of data type double: 'LB'.'

%The code below solves two equations. The initial part of the code generates 10000 values of random data U,V,x_0,y_0 which are seed values for fsolve



n=100;
%Generating 10000 values of Gaussian distributed random delay in the range
%5*10^-3 and 6*10^-3
a=5*10^-3;
b=6*10^-3;
X = (b-a).*randn(1,n) + a;

X1 = X;
d=10*10^-3;

%Generating T_1_i_A(:,1)
a=0.005;
b=5;
x = (b-a).*rand(1,n)+a;

T_1_i_A=zeros(n,30);
T_1_i_A(:,1)=transpose(x);

%Generating T_1_i_A for N=2:30 where N is no. of rounds of message exchange
for j=2:30
    T_1_i_A(:,j)=T_1_i_A(:,j-1)+X(j);
end

%Generating T_2_i_B and T_2_i_P from T_1_i_A,d,X
T_2_i_B=zeros(n,30);

for j=1:30
 T_2_i_B(:,j)=T_1_i_A(:,j)+20*10^-6+1.01*(T_1_i_A(:,j)-T_1_i_A(:,1))+15*10^-3+transpose(X);
end

T_2_i_P=zeros(n,30);

for j=1:30
 T_2_i_P(:,j)=T_1_i_A(:,j)+30*10^-6+1.105*(T_1_i_A(:,j)-T_1_i_A(:,1))+5*10^-3+transpose(X1);
end

%finding U and V
U=T_2_i_P-T_2_i_B;
V=T_1_i_A-T_1_i_A(:,1);

fun = @root2d;

%Generating 10000 random values of pre-defined clock offset in the range
%5*10^-6 and 30*10^-6
a=5*10^-6;
b=30*10^-6;
n=100;
x_0 = (b-a).*rand(1,n) + a;

%Generating 10000 random values of pre-defined clock skew in the range
%1.01 and 1.20
a1=1.01;
b1=1.20;
y_0=(b1-a1).*rand(1,n)+a1;


x0=zeros(2,n);
x=zeros(2,n);
for i=1:n
    x0(:,i)=[x_0(i); y_0(i)];
end

options = optimoptions(@fsolve, 'display', 'none');
maxN = 30;
n=100;

x = zeros(size(x0,1), n, maxN);

for N=1:maxN
    for i=1:n
        x(:,i,N)= fsolve(@(x)root2d(x, U(i,N),V(i,N),N), x0(:,i), options);
    end
    mse_theta_offset_hat_final(N)=mse(x(1,:,N),x_0);
    mse_theta_skew_hat_final(N)=mse(x(2,:,N),y_0);
end



     

   

%Plots
k=1:30;
subplot(2,1,1),plot(k,mse_theta_offset_hat_final,'r+-')
title('clock offset')
xlabel('No. of Observations')
ylabel('MSE')
grid on
subplot(2,1,2),plot(k,mse_theta_skew_hat_final,'-*')
title('clock skew')
xlabel('No. of Observations')
ylabel('MSE')
grid on
function f = root2d(x, U, V, N)
k=0;
for j=1:N
    k=k+1.5./(U-x(1)-x(2)*V+10^-4)-0.05./(U-x(1)-x(2)*V+10^-4).^2+0.1*j;
end
f(1)=k;
k=0;
for i=1:N
    k=k+(1.5*V)./(U-x(1)-x(2)*+10^-4)-0.05*V./(U-x(1)-x(2)*V+10^-4).^2+0.1*V;
end
f(2)=k;
return 
end