How to specify limits for lsqnonlin

Question

Stine Larsen 2014-9-18

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/155262-how-to-specify-limits-for-lsqnonlin

评论： Matt J 2014-9-19

I would like to specify that variables in my lsqnonlin fit are real. The other vectors in my problem are complex so I cannot use a fully real solver.

My function has the following form. r and theta are real-valued coordinate matrices. Intensity is a complex valued matrix. p is the real-valued and positive vector (at least that is what I want to specify). Hej is the function to be minimized.

Hej=  @(p)(p(1)*besselj(1,p(2)*r)+(p(4)*besselj(1,p(6)*r).*cos(theta+p(8))+p(9)*besselj(1,p(6)*r).*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1) ...
                            + (p(1)*besselj(1, p(2)*a)/besselk(1, p(3)*a)*besselk(1,p(3)*r)+...
                               (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a).*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a).*sin(theta+p(8))).*exp(1i*p(5))).*(r/a>1)...
                            -Intensity;
opts = optimoptions(@lsqnonlin,'DiffMaxChange', 0.1,'FinDiffType', 'central', 'Display','off','MaxFunEvals',2E7,'TolFun',1E-18,'TolX',1E-24,'MaxIter',4E3);
    x0 = st; % arbitrary initial guess
    lb = 0.0*ones(size(st));
    [p_estimated,resnorm,residuals,exitflag,output] = lsqnonlin(Hej,x0, lb,[], opts);

The initial guess is given as real-valued and positive vector.

So how do I specify that the only valid solution are positive and realvalued?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2014-9-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/155262-how-to-specify-limits-for-lsqnonlin#answer_152115

编辑：Matt J 2014-9-18

在 MATLAB Online 中打开

Implement the objective function as follows, splitting Hej into real and imaginary parts,

 function out=objective(p,Intensity,r,theta)
    Hej = (p(1)*besselj(1,p(2)*r)+(p(4)*besselj(1,p(6)*r).*cos(theta+p(8))+p(9)*besselj(1,p(6)*r).*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1)  + (p(1)*besselj(1, p(2)*a)/besselk(1, p(3)*a)*besselk(1,p(3)*r)+(p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a).*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a).*sin(theta+p(8))).*exp(1i*p(5))).*(r/a>1)-Intensity;
             out=[real(Hej); imag(Hej)]; %Split into real and imaginary parts
 end

and the minimization as

lsqnonlin(@(p)objective(p,Intensity,r,theta) , x0, lb,[], opts);

8 个评论
显示 6更早的评论隐藏 6更早的评论

Stine Larsen 2014-9-19

在 MATLAB Online 中打开

I managed to find a way around the first problem reformulating the objective function to:

function out=objective(p, Intensity, theta, r)

       Hej=  @(p) (p(1)*besselj(0,p(2)*r)+besselj(1,p(6)*r).*(p(4)*cos(theta+p(8))+p(9)*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1) ...
                              + (p(1)*besselj(0, p(2)*a)/besselk(0, p(3)*a)*besselk(0,p(3)*r)+...
                                 (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*sin(theta+p(8))).*besselk(1, p(7)*r).*exp(1i*p(5))).*(r/a>1)...
                                 -Intensity;
                             Hej1= @(p) real((p(1)*besselj(0,p(2)*r)+besselj(1,p(6)*r).*(p(4)*cos(theta+p(8))+p(9)*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1) ...
                              + (p(1)*besselj(0, p(2)*a)/besselk(0, p(3)*a)*besselk(0,p(3)*r)+...
                                 (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*sin(theta+p(8))).*besselk(1, p(7)*r).*exp(1i*p(5))).*(r/a>1)...
                                 -Intensity);
                             Hej2= @(p) imag((p(1)*besselj(0,p(2)*r)+besselj(1,p(6)*r).*(p(4)*cos(theta+p(8))+p(9)*sin(theta+p(8)))*exp(1i*p(5))).*(r/a<=1) ...
                              + (p(1)*besselj(0, p(2)*a)/besselk(0, p(3)*a)*besselk(0,p(3)*r)+...
                                 (p(4)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*cos(theta+p(8))+p(9)*besselj(1, p(6)*a)/besselk(1, p(7)*a)*sin(theta+p(8))).*besselk(1, p(7)*r).*exp(1i*p(5))).*(r/a>1)...
                                 -Intensity);
                             out= @(p) [Hej1; Hej2]; %Split into real and imaginary parts

But now I run into a second problem, lsqnonlin only accepts inputs in double format and I get the following error:

Error using lsqnonlin (line 233)
LSQNONLIN requires all values returned by user functions to be of data type double.
Error in Mode_decomposition_2modes_field_real_imag (line 191)
    [p_estimated,resnorm,residuals,exitflag,output] = lsqnonlin(@(p) objective(p,Intensity,r,theta), x0, lb,[], opts);

Stine Larsen 2014-9-19

Using a cell array in 'out' this is solved.

Matt J 2014-9-19

Way, way, way over-complicated. You have functions within functions within cells, whereas my original proposal was just two lines of numeric MATLAB operations...

请先登录，再进行评论。

Answer 2

Roger Wohlwend 2014-9-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/155262-how-to-specify-limits-for-lsqnonlin#answer_152086

I am afraid you cannot order the optimizing function to search for a real solution. If you want a real solution you have to make sure that your function Hej returns only real values. Since you use lsqnonlin as optimizing function your goal is to minimize the sum of the squared elements of the vector Hej. In my opinion, it does not make sense to minimize the sum of squared complex numbers. It would make more sense to minimize the sum of the absolute values of the complex numbers. So perhaps you should adjust Hej in a way that it returns the absolute values of the complex numbers instead of the complex numbers. If you do that, Hej returns only real values. As a consequence lsqnonlin should search for real solutions only.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Stine Larsen 2014-9-18

Thank you for the suggestion. I will try it out and see if it gives a physical valid solution. I would still like to consider the full complex structure in case anyone has other suggestions.

请先登录，再进行评论。