Use optimisation toolbox with external programme

Hi,

I want to use the optimisation toolbox, in particular the lsqcurvefit function, to fit the response of an external C programme to experimental data. For testing purposes I have written a C programme which computes f(x) = p*sqr(x), where p is my parameter to be optimised. However, the optimsation alogorithm does not perform any optimsation of my parameter p. It stops after two function evaluation:

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point 
is less than the default value of the function tolerance.

I checked the first-order optimality and it is zero:

firstorderopt: 0

This is my function to the external programme:

% ---------------  external function (m-file) ---------------------
function y = myexternal(p,xdata)   % p = parameter to be optimized, xdata = measuring points

% write parameter to file
writeToFile(p)

% run my external
system('myprog');

% read result from file
data = textscan(resultsfile,'%f');
y=data{1};

end
% ---------------------- end external function (m-file) ----------------------------

I run my optimisation (assuming xdata and ydata is definied) via:

p0 = [0.3];
[p,resnorm,~,exitflag,output,lambda] = lsqcurvefit(@myexternal,p0,xdata,ydata)

Interesting to note that if I replace in the function myexternal() the line "y=data{1};" by "y=p*xdata.*xdata" everything works fine. I also compared the numbers in the computed array y in both cases and it contains the same values.

Do you have any help what might be the mistake? Help is highly appreciated: :-)

Maik

Answer 1

[Maik]

Dear Matt,

I got it. The reason why it didn't work was related to the output accuracy of the matlab script and the external C programme.

This is the output for the matlab script (using 16 decimal places for double variables)

% writing paramater as output
fprintf(fp_parameter_file,'%.16g\n',p);

% reading result of C programme
data = textscan(fp_results_file,'%f %f');

And for the C programme

 % reading parameter input of matlab script
 fscanf(fp_parameter_file, "%lf\n",&p);

 % writing results file
 fprintf(fp_results_file, "%.16g\n", result);

Thanks for your support!

Answer 2

[Matt J]

This isn't quite what you asked, but if your model function has a linear form in p

   F(p,xdata)=p*sqr(xdata)

and you have no bounds on p, then it makes no sense to be using lsqcurvefit. The result has the trivial analytical solution

    p= reshape(sqr(xdata),[],1)\ydata(:);

As for the code you've shown, it's hard to interpret, because you haven't explained anywhere which part is computing sqr(x). However, one obvious problem is that you are not using xdata anywhere inside the version of myexternal() that you've posted.

[Maik]

Hey Matt,

thanks for your reply. The example provided is just a very simple test case. Basically, I want to couple the optimisation toolbox of Matlab with a finite-element (FE) programme. I want to optimise the material parameters of numerical model to best fit experimental data. I was also wondering if the unused xdata might be the problem. But than I found an older post ( link ) where someone used the optimisation toolbox together with an FE programme. However, is question was addressing the multistart feature. So it is working, but I did not get further explanations on how to do it.

Regarding your question: The C programme myprog reads a file where the parameter p is stored, computes an array of function values for an given intervall and writes to a file, which will be read by myexternal() in stored in the array y.

In summary, myexternal() takes the parameter p returns and an array of values. This it is basically the same as defining "y=p*xdata.*xdata" instead of "y=data{1}". At least for understanding, but obviously it is not. :-(

Hopefully this makes things clearer.

[Matt J]

In summary, myexternal() takes the parameter p returns and an array of values. This it is basically the same as defining "y=p*xdata.*xdata" instead of "y=data{1}".

Well you should run myexternal() in isolation and see if it is indeed returning the equivalent of p*xdata.*xdata.

[Maik]

I just checked it. It really returns the same value. I'm running out of ideas at the moment. :-(

[Matt J]

Add this to the last line of your objective function

discrepancy = norm(y-p*xdata.^2)/norm(p*xdata.^2)

Run the solver and see if it always generates small discrepancies every time the function is called.

[Maik]

Ich justed check the discrepancy and it is 1.0610e-16. I would conclude that the external programme returns the same result. Obviously, Matlab is performing some steps which are not covered by the external-programme version. :-/

I just checked the number of function calls. Its two times at each iteration , which i assume is ok in order to compute the gradient numerically.

[Maik]

However, I just realised that the if I output the parameter p inside the function to be optimised, at each iteration step, I only get the newly computed parameter, but I do not get the pertuabted parameter, which is needed for the finite difference method.

[Matt J]

If everything you say is true, the only thing I can think of is that your initial p0 is already optimal, and therefore lsqcurvefit stops immediately. Either that, or else sqr(x)=0, which would imply that all values of p are optimal.

Since your practice problem only has a single parameter, p, you could plot the objective function

f=@(p) norm(y-myexternal(p,xdata))^2

as a function of p to make sure it varies non-trivially around p0 (it should be a simple quadratic parabola). I would recommend posting such a plot here so that we have a feel for what's happening, but make sure you generate the plot for a case when the fit is failing.

You could also try minimizing the above with other solvers, like fminsearch, just to see if the results are different.