multistart using lsqcurvefit or lsqnonlin -凯发k8网页登录

step 1. create the objective function.

write an anonymous function that takes a data matrix xdata with n rows and two columns, and returns a response vector with n rows. the function also takes a coefficient matrix p, corresponding to the coefficient vector (a,b,c,d).

fitfcn = @(p,xdata)p(1)   p(2)*xdata(:,1).*sin(p(3)*xdata(:,2) p(4));

step 2. create the training data.

create 200 data points and responses. use the values a=-3,b=1/4,c=1/2,d=1. include random noise in the response.

rng default % for reproducibility
n = 200; % number of data points
preal = [-3,1/4,1/2,1]; % real coefficients
xdata = 5*rand(n,2); % data points
ydata = fitfcn(preal,xdata)   0.1*randn(n,1); % response data with noise

step 3. set bounds and initial point.

set bounds for lsqcurvefit. there is no reason for d to exceed π in absolute value, because the sine function takes values in its full range over any interval of width 2π. assume that the coefficient c must be smaller than 20 in absolute value, because allowing a high frequency can cause unstable responses or inaccurate convergence.

lb = [-inf,-inf,-20,-pi];
ub = [inf,inf,20,pi];

set the initial point arbitrarily to (5,5,5,0).

p0 = 5*ones(1,4); % arbitrary initial point
p0(4) = 0; % ensure the initial point satisfies the bounds

step 4. find the best local fit.

fit the parameters to the data, starting at p0.

[xfitted,errorfitted] = lsqcurvefit(fitfcn,p0,xdata,ydata,lb,ub)

local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

xfitted = 1×4
   -2.6149   -0.0238    6.0191   -1.6998

errorfitted = 28.2524

lsqcurvefit finds a local solution that is not particularly close to the model parameter values (–3,1/4,1/2,1).

step 5. set up the problem for multistart.

create a problem structure so multistart can solve the same problem.

problem = createoptimproblem('lsqcurvefit','x0',p0,'objective',fitfcn,...
    'lb',lb,'ub',ub,'xdata',xdata,'ydata',ydata);

step 6. find a global solution.

solve the fitting problem using multistart with 50 iterations. plot the smallest error as the number of multistart iterations.

ms = multistart('plotfcns',@gsplotbestf);
[xmulti,errormulti] = run(ms,problem,50)

multistart completed the runs from all start points.
all 50 local solver runs converged with a positive local solver exit flag.

xmulti = 1×4
   -2.9852   -0.2472   -0.4968   -1.0438

errormulti = 1.6464

multistart finds a global solution near the parameter values (–3,–1/4,–1/2,–1). (this is equivalent to a solution near preal = (–3,1/4,1/2,1), because changing the sign of all the coefficients except the first gives the same numerical values of fitfcn.) the norm of the residual error decreases from about 28 to about 1.6, a decrease of more than a factor of 10.

formulate problem for lsqnonlin

for an alternative approach, use lsqnonlin as the fitting function. in this case, use the difference between predicted values and actual data values as the objective function.

fitfcn2 = @(p)fitfcn(p,xdata)-ydata;
[xlsqnonlin,errorlsqnonlin] = lsqnonlin(fitfcn2,p0,lb,ub)

local minimum possible.
lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

xlsqnonlin = 1×4
   -2.6149   -0.0238    6.0191   -1.6998

errorlsqnonlin = 28.2524

starting from the same initial point p0, lsqnonlin finds the same relatively poor solution as lsqcurvefit.

run multistart using lsqnonlin as the local solver.

problem2 = createoptimproblem('lsqnonlin','x0',p0,'objective',fitfcn2,...
    'lb',lb,'ub',ub');
[xmultinonlin,errormultinonlin] = run(ms,problem2,50)

multistart completed the runs from all start points.
all 50 local solver runs converged with a positive local solver exit flag.

xmultinonlin = 1×4
   -2.9852   -0.2472   -0.4968   -1.0438

errormultinonlin = 1.6464

again, multistart finds a much better solution than the local solver alone.