MATLAB Answers

Help with global optimization problem

6 views (last 30 days)
adullak on 16 Sep 2021
Commented: Star Strider on 17 Sep 2021
So I'm trying to fit a model to data using global search (I also tried all of the global optimization solvers without a better result). Basically I am trying to fit a signal from an instrument. This signal is supposed to be composed of three overlapped peaks, and I am trying to deconvolve the signal into those 3 individual peaks so I can analyze them further one by one. The way I'm trying to do this is by fitting the signal to a sum of three functions, so called Fraser Suzuki functions (which I read somewhere are the same as the lognormal function), each of which has 4 parameters, so 12 parameters overall.
I created a toy version of my problem, in which instead of using real data, mock data is generated using said sum of these functions. I then try to use global optimization algorithms to minimize the sum of square differences to get a result and then compare it to the parameters found by the minimizations to the ones used to generate the data.
The problem is that the parameters determined by minimization are not the same as the original ones, and in theory one should be able to determine the original parameters by minimization, since the data were created with the same functions I am trying to fit to them. I got to get good results by imposing really stringent upper and lower bounds, but I have to analyze a lot of signals and trying to make this as automatic as possible is the reason I'm using code in the first place.
I'd be glad to read some of your input.
Thank you very much in advance.
PS: I ended up adding the pattern search solver since it returned a way better result that globalsearch.
Here's the code
%Parameters to create data
pa = [0.0363 0.0682 560.0000 55.1830 0.0597 -0.2543 632.3738 35.0086 0.0118 0.9939 609.8258 160.0000];
%independent variable vector
Temp = (400:1000)';
%create data
fs_data = fsm(Temp,pa);
%set lower and upper bounds
lb = [0, -.5, 0, 0,...
0, -.5, 0, 0,...
0, -.5, 0, 0];
ub = [0.1, 0.5, 800, 200,...
0.1, 0.5, 800, 200,...
0.1, 0.5, 800, 400];
%Set initial guess
p0 = [0.0363 0.0682 560.0000 55.1830 0.0597 -0.2543 632.3738 35.0086 0.0118 0.9939 609.8258 160.0000];
%Set objective function
obj_fun =@(parameters) sum((fsm(Temp,parameters)-fs_data).^2);
%Set optmization problem
rng default
[fsfit,fval] = patternsearch(obj_fun,p0,[],[],[],[],lb,ub)
% % Set optmization problem
% rng default % For reproducibility
% opts = optimoptions(@fmincon,'Algorithm','sqp');
% gs = GlobalSearch;
% problem = createOptimProblem('fmincon','x0',p0,'objective',obj_fun,'lb',lb,'ub',ub,'options',opts);
% % fit model to data
% [fsfit,fval] = run(gs,problem)
difference = pa - fsfit
hold on
hold off
%Define Fraser Suzuki Function
function y = fs(T,p)
in = 2 * p(2) * ((T-p(3)) / p(4));
out = - log(2) / p(2)^2;
y = p(1) * exp(out .* (log(1+in)).^2);
%Define Fraser Suzuki mixture
function y = fsm(T,pm)
y = fs(T,pm(1:4)) + fs(T,pm(5:8)) + fs(T,pm(9:12));

Answers (1)

Star Strider
Star Strider on 16 Sep 2021
The Golbal Optimization Toolbox functions (as much as I like them) may not be appropriate for this. See if the rica function (that performs an independent components analysis decomposition) will separate them. If so, it will then be possible to estimate the parameters describing the individual signals.
Star Strider
Star Strider on 17 Sep 2021
My pleasure!
Without your actual (or simulated) data, a more precise response is not possible.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by