MultiStart with lsqnonlin, Problem-Based

This example shows how to fit a function to data using lsqnonlin together with MultiStart in the problem-based approach.

Many fitting problems have multiple local solutions. MultiStart can help find the global solution, meaning the best fit.

The model is

$$y=a+b{x}_1\sin (c{x}_2+d)$$,

where the input data is $$x=(x_1,x_2)$$, and the parameters a, b, c, and d are the unknown model coefficients.

Create Problem Data

Most problems involve data from measurements. For this problem, create artificial data including noise. Create 200 data points and responses. Specify the values $$a=-3$$, $$b=1/4$$, $$c=1/2$$, and $$d=1$$.

rng default % For reproducibility
fitfcn = @(p,xdata)p(1) + p(2)*xdata(:,1).*sin(p(3)*xdata(:,2)+p(4));
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

Create Optimization Variables and Problem

The optimization variables are the model coefficients. Set bounds for the variables as you create them. The variable $$d$$ does not need to exceed $$\pi$$ in absolute value, because the sine function takes values in its full range over any interval of width $$2\pi$$. 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.

a = optimvar("a");
b = optimvar("b");
c = optimvar("c",LowerBound=-20,UpperBound=20);
d = optimvar("d",LowerBound=-pi,UpperBound=pi);
prob = optimproblem;

Create Objective Function

Create the objective function as the sum of squared differences between the model response and the data.

resp = a + b*xdata(:,1).*sin(c*xdata(:,2) + d);
prob.Objective = sum((resp - ydata).^2);

Create Initial Point and Solve Problem

The initial point is a structure with values for the coefficients. Arbitrarily set the initial point to [5,5,5,0].

x0.a = 5;
x0.b = 5;
x0.c = 5;
x0.d = 0;
[sol,fval] = solve(prob,x0)

Solving problem using lsqnonlin.

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.

sol = struct with fields:
    a: -2.6149
    b: -0.0238
    c: 6.0191
    d: -1.6998

fval = 
28.2524

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

Search for Improved Solution Using MultiStart

To search for a better solution, create a MultiStart object. Plot the best function value as the solver iterates.

ms = MultiStart(PlotFcns=@gsplotbestf);

Call solve again, this time using ms. Specify 50 start points for MultiStart.

rng default % For reproducibility
[sol2,fval2] = solve(prob,x0,ms,MinNumStartPoints=50)

Solving problem using MultiStart.

MultiStart completed the runs from all start points. 

All 50 local solver runs converged with a positive local solver exitflag.

sol2 = struct with fields:
    a: -2.9852
    b: -0.2472
    c: -0.4968
    d: -1.0438

fval2 = 
1.6464

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