Solving nonlinear equation using lsqnonlin

Hi

I am trying to solve a set of non-linear equations using lsqnonlin. Initially I was using 'levenberg-marquardt', but I was getting arbitrary results than what I was supposed to get. Then I tried solving using 'lsqnonlin' by defining some upper and lower bounds for the variables. However, the variables are now clinging to either the upper or the lower bound value. Therefore, the reults once again are not correct. What might be the reason for such behaviour?

Here is the code that I am using: (There are 36 to 56 variables. I am showing just 2 for example)

lb=[0.5,-0.5];
ub=[1.5,0.5];
x0=[1,0];
options=optimoptions(@lsqnonlin,'Algorithm','trust-region-reflective'); 
[x,res]=lsqnonlin(fun,x0,lb,ub,options);

[Ameer Hamza]

This is likely caused by your objective function 'fun'. It may be unbounded, so the optimal solution is clinging to the extreme limits.

[Anuj Kumar Sahoo]

Hi Ameer

Thanks for your response. I chnaged the objective function as per your suggestion and the variables are no longer clinging to the upper and lower bounds.

However, there is another problem I hope you can help. I have a system which I am trying to model. When the system is ideal, I know the values of the variables. When I use those values in my objective function I get good results. However, when I am trying to findout the values of those unknown coefficients/variables of the same ideal system using 'lsqnonlin' they are coming far from the expected results. That is why I tried setting limits to the variables.

After implementing your suggestion the coefficients are no longer at their extreme limits, but still they are not giving exact result. Shouldn't I get the coefficients very close or even same as the ideal system coefficients, since I am setting same values as my initial guess? What might be the reason for this? I hope I have made myself clear.

Thanks in advance.

Answer 1

[Walter Roberson]

When I use those values in my objective function I get good results. However, when I am trying to findout the values of those unknown coefficients/variables of the same ideal system using 'lsqnonlin' they are coming far from the expected results.

There is no possible algorithm that can find the global minimum of arbitrary "black box" functions (a "black box" function is one that has to just be executed but cannot be effectively examined.) Literally not possible -- it has been proven mathematically.

Different minimizer algorithms have different kinds of functions that they do a good job on. And, conversely, different kinds of situations that they tend to get stuck in.

lsqnonlin often does much better than I expect, but there are some situations it does not do well on.

One example of a situation lsqnonlin does not do well on, is the sum of two gaussians that are slightly different formulas, and which are not well separated. In such a situation, lsqnonlin tends to drive one of the two gaussians to effectively take over all of the fitting, with the other one being driven indefinitely wide -- effectively the second one gets driven to be a constant term. If you start too close to one side the first gaussian takes over; if you start too close to the other side, the other gaussian takes over. You might need to start quite close to the true positions for lsqnonlin to find the correct values. lsqnonlin is not able to find the correct positions itself because every time it moves away from the "one gaussian plus effectively constant" situation, the fit seems to get much much worse.

Another case that optimizers have trouble with is a central region that is narrow and deep, with asymptoptic descent on either side of that: if the optimizer happens to land outside the edge of the central region, then the slope is always decreasing away from the central region, so the optimizers tend to chase that out until infinity or the bounds. The optimizer would have to deliberately "go uphill" hoping that there might just happen to be a well to be discovered.

These kinds of problems are often difficult to deal with unless you can find an optimizer that is designed for the particular shape of function you are dealing with. There are algorithms for sum of gaussians, for example (but I have never researched to find out how they work.)

Answer 2

[Ameer Hamza]

As Walter already explained that there is no guaranteed way to get a globally optimum solution for an arbitrary problem using a numerical method. Because of the way they are formulated, gradient-based methods can at best reach a locally optimal solution depending on the initial guess. Some metaheuristic optimizers, such as genetic algorithm ga(), or particle swarm optimizer particleswarm() have a higher chance of reaching a global solution, but still, the globally optimum solution can never be guaranteed. However, you can increase the probability in several ways. Global Optimization toolbox gives necessary tools for that.

For example, see GlobalSearch(): https://www.mathworks.com/help/gads/globalsearch.html or multistart(): https://www.mathworks.com/help/gads/multistart.html functions. They provide a systematic way to run the common optimizers, such as fmincon(), with several different starting points hoping that this will lead to a global solution. Similarly, check ga(): https://www.mathworks.com/help/gads/ga.html and particleswarm(): https://www.mathworks.com/help/gads/particleswarm.html.

[Anuj Kumar Sahoo]

Hi Ameer and Walter

Thanks a lot to both of you for your suggestions. It helped me to look into my problem more deeply. In order to findout the cause of error in my result I took a small portion of my system, simulated it and tried to fit the system using some simple equations. To verify my equations I changed some parameters from their ideal value in my simulaion set up and tried to find out whether I can retrieve those modified parameters from the equations or not.

This is part of the code which I am using:

%%% %%%%%%
p1=x1*s21*exp(-1i*x3);
p2=x2*s31*exp(-1i*x4);
pin1=(abs(p1*s21+p2*s31))^2;
pin2=(abs(p2*s21+p1*s31))^2;
pin3=(abs(p2*p1+p1*p2))^2;
pin4=(abs(p1*p1+p2*p2))^2;
F(1)=pin11a-a;
F(2)=pin12a-b; 
F(3)=pin13a-c; 
F(4)=pin14a-d; 
%%%%%%%%%

s21 and s31 are the ideal circuit parameters (complex valued). x1, x2, x3, x4 are the variables that I have introduced to model their amplitude and phase imperfections. a, b, c, d are the readings I obtained from simulator. So, my idea is if modeled correctly F(1), F(2), F(3), F(4) must be zero/minimum.

However, after solving this I found out that the variables asoociated with phase imperfections (x3, x4) are not changing at all and returning the initial guess value. Therefore, the variables x1, x2 are adjusting themselves for best fit and turning out be incorrect. Since, this is an important logic of my original code I believe this may be causing some error.

I have tried to explain my problem in detail. I hope it is clear to you. Kindly help why the phase variables are not getting modeled.

Thanks in advance

[Ameer Hamza]

Without being able to run this example, it is difficult to say anything. Can you attach a small example, which we can run. It will make it easier to suggest anything.

[Anuj Kumar Sahoo]

Hi Ameer

Thanks for your interest. I am attaching the code herin.

%%%%% main function %%%%%
G=[0.0002 0.9279 0.0015 0.8832]; % readings from simulator
fun=@(x)BLC_model_call(x,G); 
%%%%%%%initial guess %%%%%%%%%%%
x0(1)=1;
x0(2)=1;
x0(3)=-0.1;
x0(4)=-0.1;
opts.Algorithm='levenberg-marquardt'; 
opt.Jacobian='romberg';
[x,res1,residual]=fsolve(fun,x0,opts);
%%%%%%%%%%% end %%%%%%%
%%%%% objective function %%%%%
function F=BLC_model_call(x,G)
x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=G(1);
b=G(2);
c=G(3);
d=G(4);
%%%%%% derivation of ideal system parameters %%%%%%%%%%
z0=50;z1=50;z2=35.35;
freq=2e9;
fr=2e9;
th1=1i*pi*freq/(4*fr);
th2=1i*pi*freq/(2*fr);
blce=[1 0;1/(z1*coth(th1)) 1]*[cosh(th2) z2*sinh(th2);(1/z2)*sinh(th2) cosh(th2)]*[1 0;1/(z1*coth(th1)) 1]; % abcd blc even
blco=[1 0;1/(z1*tanh(th1)) 1]*[cosh(th2) z2*sinh(th2);(1/z2)*sinh(th2) cosh(th2)]*[1 0;1/(z1*tanh(th1)) 1]; % abcd blc odd
%% even gamma and transmission
gblce=(blce(1)+blce(3)/z0-blce(2)*z0-blce(4))/(blce(1)+blce(3)/z0+blce(2)*z0+blce(4));
tblce=2/(blce(1)+blce(3)/z0+blce(2)*z0+blce(4));
%% odd gamma and transmission
gblco=(blco(1)+blco(3)/z0-blco(2)*z0-blco(4))/(blco(1)+blco(3)/z0+blco(2)*z0+blco(4));
tblco=2/(blco(1)+blco(3)/z0+blco(2)*z0+blco(4));
%% blc s-parameters
blcs21=0.5*(tblce+tblco);
blcs31=0.5*(tblce-tblco);
blcs24=blcs31;
blcs34=blcs21;
%%%%% modeling system parameters %%%%%%
pina1=x1*blcs21*exp(-1i*x3);
pina2=x2*blcs31*exp(-1i*x4);
%%%% deriving power %%%%%%
pin11a=(abs((pina1*blcs21*1+pina2*blcs24*1)))^2;
pin12a=(abs((pina2*blcs34*1+pina1*blcs31*1)))^2;
pin13a=(abs((pina1*pina1*1+pina2*pina2*1)))^2;
pin14a=(abs((pina2*pina1*1+pina1*pina2*1)))^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F(1)=pin11a-a; 
F(2)=pin12a-b; 
F(3)=pin13a-c; 
F(4)=pin14a-d; 
%%%%%% end %%%%%

The result I am getting:

x(1)=0.968435179363055
x(2)=0.967701516736763
x(3)=-0.0999996983986709
x(4)=-0.100000404052192

The result I should get:

x(1)=0.97772
x(2)=0.95063
x(3)=-0.3839
x(4)=-0.3752

I hope this code is sufficient for you to analyze. I hope you can figure out some solution for this problem.

Thanks a lot for your time and help.

[Walter Roberson]

You are using abs() . However, the derivative of abs() is discontinuous, which is a problem for fsolve(). Especially with you using abs()^2 you might be able to rewrite the calculations without using abs() -- in particular, abs(a)^2 == a.*conj(a)

[Anuj Kumar Sahoo]

Hi Walter

Thanks a lot for your time.That was an insightful observation from you. I changed my code as per your suggestion and represented "abs(a)^2" in the form of "a.*conj(a)" . Unfortunately the problem still persists and the variable with phase imperfections are not getting modelled. They are returning the initial guess which I am using.

Here is the modified section in my code:

pin11a=((pina1*blcs21+pina2*blcs24)*conj(pina1*blcs21+pina2*blcs24));

pin12a=((pina2*blcs34+pina1*blcs31)*conj(pina2*blcs34+pina1*blcs31));

pin13a=((pina1*pina1+pina2*pina2)*conj(pina1*pina1+pina2*pina2));

pin14a=((pina2*pina1+pina1*pina2)*conj(pina2*pina1+pina1*pina2));

Is it something to do with the presentation of pina1 as,

pina1=x1*blcs21*exp(-1i*x3);

I have also tried presenting this way,

pina1=x1*blcs21*(cosd(x3)-1i*sind(x3));% argument in terms of degree

But there was no change in the result.

Kindly help.

Thanks.

[Anuj Kumar Sahoo]

Hi Walter

I represented pina1 and pina2 as

pina1=blcs21*(x1+1i*x3);

pina2=blcs31*(x2+1i*x4);

Now I am getting the desired result.

Thanks for helping me.

[Walter Roberson]

The editor points out that you assign to gblce and gblco but you do not use those variables. But on the line after you assign to gblce you use blce and on the line after you assign to gblco you use blco . Is it possible that you should have been using the g* version of the variables on those lines?

[Anuj Kumar Sahoo]

Hi Walter

I am sorry for such confusing variable names. Actually, these are standard equations to derive the system parameters. The parameters gblce and gblco are used for calcualting reflected signal. Similary tblce and tblco are used for transmitted signal calculation. In my case I only need transmitted signal, so I have used tblce and tblco. So that part is correct.

However, I have found out that the main problem in my code might be (abs(expression))^2. Which fsolve is not able to solve. So, I will be grateful if you could suggest any alternate command or code to resolve this limitation.

FYI: I have tried this, abs(a)^2 == a.*conj(a). But it didn't work.

Thanks