How can I use fmisearch() to evaluate the parameters in a matrix?

I have the following equation:

f(t)=integral(g(t,x)*h(x) dx)

I discretized the problems because the symbolic expression can't be evaluate (too complicated), so:

where:

f(t) is a vector length=n;

g(t,x) is a matrix n rows and m columns;

h(x) should be a vector length=m. However it depends on 3 parameters.

I want to determine those 3 parameters such that f(t)-integral()=0.

If the 3 parameters are known I use trapz or Simpson method to evaluate the integral.

Do you know a method to apply fminsearch() to evaluate the 3 parameters?

Thank you in advance

Answer 1

[Star Strider]

If the objective is to determine the parameter at a zero-crossing (root), use fsolve or fzero, not fminsearch.

Try something like this —

g = @(t,x) sin(2*pi*x.*t);                          % Create Function
h = @(x) exp(x);                                    % Create Function

a = -5;                                             % Define Lower Integration Limit
b = +5;                                             % Define Upper Integration Limit
f = @(t) integral(@(x)g(t,x).*h(x), a, b)           % Integral Funciton Call (Objective Fucntion)

f = function_handle with value:
    @(t)integral(@(x)g(t,x).*h(x),a,b)

[T,Fval] = fsolve(f, 42)                            % Estimate Parameter

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

T = 143.2500

Fval = -1.3207e-10

.

[3Nz0]

thanks for your answer.

The problem is that g(t,x) is a a function more complicated:

x(t,mw)=t./(M_star*k_tdd*(mw.^2));

gx(t,mw)=-x+(x.^0.5).*(x+((pi*x).^0.5)+pi).^0.5;

U(t,mw)=((1/k_tdd)*(1./mw.^3).*t)+(M_star*gx./mw);

F_sr(t,mw) = (8/pi^2)*symsum(1/((2*I-1)^2)*exp(-(2*I-1)^2 .*U), I, 1, inf);

w(mw) = ((mw./M0).^(a+1)).*exp(-(mw/M0).^b);

this is the list of all functions:

I have to solve the integral of F*w. The parameters I want to find are: a, b, M0 in w(). They are such that the result (a vector) is equal to another know vector (experimental data).

For this reason I discretized the problem.

[Star Strider]

As always, my pleasure!

I have no idea what you want to do with this, and several parameters are missing, as well as a syms call initially.

I also do not know what the arguments are to integral or what ‘f(t)’ is in:

‘I want to determine those 3 parameters such that f(t)-integral()=0’

(There is also a symsum call, however I do not understand what needs to be summed.) This approach does everything numerically, since it is likely that this is more efficient than a symbolic solution for this problem. It may be that this can be set up symbolically, then solved numerically. I will need more information in that regard.

Please supply those, and I will do my best to get this working.

x = @(t,mw) t./(M_star*k_tdd*(mw.^2));
gx = @(t,mw) -x(t,mw)+sqrt(x(t,mw)).*sqrt(x(t,mw)+(sqrt(pi*x(t,mw)))+pi);
U = @(t,mw) ((1/k_tdd)*(1./mw.^3).*t)+(M_star*gx(t,mw)./mw);
F_sr = @(t,mw)  (8/pi^2)*integral(@(I) 1./((2*I-1).^2).*exp(-(2*I-1).^2 .*U(t,mw)), 1, inf);
% % w = @(p,mw) = ((mw./M0).^(a+1)).*exp(-(mw/M0).^b);
w = @(p,mw) = ((mw./p(3)).^(p(1)+1)).*exp(-(mw/p(3)).^p(2));            % p(1) = aa,  p(2) = b,  p(3) = M0
f = @(p,t) integral(@(mw)G_sr(t,mw).*h(mw), a, b)                       % Integral Funciton Call (Objective Fucntion)
t = 42;
[P,Fval] = fsolve(@(p)f(p,t), rand(3,1))                            % Estimate Parameter

.

[3Nz0]

Sorry, I just copied the equations from an mfile where I tried to solve it in a symbolic way.

The results from the integral should be a curve equal to f(t).

f(t) comes from the experimental data

[Star Strider]

First, a representative vector for ‘ft’ would helpl.

Second, please explain the symsum call, and please supply the syms declarations.

I want to help you get this working, however I need this informaiton in order to do that.

[3Nz0]

The 1st part of the mfile I wrote in symbolic way is:

G0N=2*10^5;

Me=18500;

M_star=160000;

k_tdd=1.05*10^(-15);

k_rouse=2*10^(-12);

beta=2.3;

syms I t mw

x(t,mw)=t./(M_star*k_tdd*(mw.^2));

gx(t,mw)=-x+(x.^0.5).*(x+((pi*x).^0.5)+pi).^0.5;

U(t,mw)=((1/k_tdd)*(1./mw.^3).*t)+(M_star*gx./mw);

F_sr(t,mw) = (8/pi^2)*symsum(1/((2*I-1)^2)*exp(-(2*I-1)^2 .*U), I, 1, inf);

wGEX(mw) = ((mw./M0).^(a+1)).*exp(-(mw/M0).^b);

w(mw)=wGEX./mw;

F_w(t,mw)=F_sr.*w(mw);

G=int(F_w, mw, 0, inf);

the results is: a=11.349; b=3.366; M0=247300;

the experimental data (Gt):

t= 100 50 25 10 5 2.63 1.03 0.526 0.244 0.0991 0.0516 0.0251 0.00996 0.00495 0.00252 0.001

G= 124 879 6110 31900 64000 94400 131000 151000 169000 184000 193000 201000 209000 215000 235000 303000

Mw= 129508.9 137732.8 145580 157710.6 172968.9 190874.8 210634.5 233875.1 248726.4 258086.2 264520.8 269451.9 274474.9 281318.2 286562.4 293707.1 301029.9 306641.6 308535.3 314286.8 318180.7 322122.7 328127.6 330154 342577.9 351119.2 359873.4 366582 368845.9 373415.7 375721.7 380376.7 387467.5 392268 399580.5 402048.2 407029.3 409543 412072.1 414617 419753.8 424954.3 427578.7 432876.1 435549.4 446408.7 451939.4 460364.3 466067.9 468946.2 474756.1 477688.1 486592.9 501804.7 520688 550353.5 581709.1 603599.2

I made another mfile where there are no function, mw and t are vectors and I use simpson method to evaluate the integral (only if a, b , M0 are known!!!).

Sorry if I'm not clear.

Thank you again

[3Nz0]

In the attachment there is the mfile with the discretized problem.

[Star Strider]

This runs without error, however it appears to be an extremely difficult problem.

I am currently runnig it with GlobalSearch, although MulitStart would also be an option, so consider that as well, if you have Global Optimization Toolbox.

I replaced the symsum call with an integral call because the summation does not have a closed-form solution, and so prevents code created with matlabFunction from running successfully in the numerical parameter estimation.

G0N=2E5;
Me=18500; 
M_star=160000;
k_tdd=1.05E-15;
k_rouse=2E-12; 
beta=2.3; 
% syms I t mw M0 a b 
digits(5)
x = @(t,mw) t./(M_star*k_tdd*(mw.^2));
gx = @(t,mw) -x(t,mw)+sqrt(x(t,mw)).*sqrt(x(t,mw)+(sqrt(pi*x(t,mw)))+pi);
U = @(t,mw) ((1./k_tdd).*(1./(mw.^3)).*t)+(M_star.*gx(t,mw)./mw);
F_sr = @(t,mw) (8/pi^2)*integral(@(I) 1./((2*I-1).^2) .* exp(-(2*I-1).^2 .* U(t,mw)), 1, inf, 'ArrayValued',1)
wGEX = @(mw,M0,a,b) ((mw./M0).^(a+1)).*exp(-(mw/M0).^b);
w = @(mw,M0,a,b) wGEX(mw,M0,a,b)./mw;
F_w = @(t,mw,M0,a,b) F_sr(t,mw).*w(mw,M0,a,b);
G = @(t,M0,a,b) integral(@(mw)F_w(t,mw,M0,a,b), 0, inf, 'ArrayValued',1)

G_fcn = @(t,M0,a,b) G(t,M0,a,b);

 t = [100  50  25  10  5  2.63  1.03  0.526  0.244  0.0991  0.0516  0.0251  0.00996  0.00495  0.00252  0.001 ];
 
 G = [124  879  6110  31900  64000  94400  131000  151000  169000  184000  193000  201000  209000  215000  235000  303000 ];
 
 G_fcnp = @(p,t) G_fcn(t,p(1),p(2),p(3));
 
 ObjFcn = @(p) norm(G - G_fcnp(p,t))
 gs = GlobalSearch;
 problem = createOptimProblem('fmincon', 'x0',randn(3,1)*1E-3, 'objective',ObjFcn);
 P = run(gs,problem)
 
%  [P,resnorm] = lsqcurvefit(G_fcnp, randn(3,1), t, G)
 
 G_out = G_fcnp(P,t)
 
 figure
 plot(t, G, 'p')
 hold on
 plot(t, G_fcnp(P,t), '-r')
 hold off
 grid
 

I left the lsqcurvefit call in (although commented-out) so use it with the value of ‘P’ returned by GlobalSearch to get other statistics. If you have the Statistics and Machine Learning Toolbox, consider using fitnlm with ‘G_fcnp’ is also an option, again using the returned value of ‘P’ as the initial parameter estimates, in order to get statistics on the parameters and the fit.

I am stopping here for the time being, and will edit this to post whatever GlobalSearch produces as a parameter estimate vector, if it succeeds.

EDIT — (30 Apr 2021 at 20:44)

After two hours of running using:

 ms = MultiStart;
 problem = createOptimProblem('fmincon', 'x0',randn(3,1)*1E-3, 'objective',ObjFcn);
 P = run(ms,problem,50)

instead of GlobalSearch (that ran for a while without finding the global minimum), the result is —

MultiStart stopped without completing the runs from all start points.

1 out of 50 local solver runs exceeded the iteration limit (problem.options.MaxIterations) or 
the function evaluation limit (problem.options.MaxFunctionEvaluations).
None of the 50 local solver runs converged with a positive local solver exit flag.

P =

       1.7273
     -0.10245
      0.10054

I will give it a greater range of parameters, run it for a while, and see if it converges (although I am not confident that it will).

[3Nz0]

Wow, thank you so much for your help.

In the function definitions there are some errors:

wGEX = @(mw,M0,a,b) 1.9*(((mw./M0).^(a+1)).*exp(-(mw/M0).^b));

G_fcn = @(t,M0,a,b) G0N.*(G(t,M0,a,b)).^beta;

Further:

-I think is better to use P0=[10 3 200000] (initial values near to the final solution).

-maybe, the definition of ObjFcn shoul be chi_squared minimization: ObjFcn = @(p) (1./n_exp).*sum(((G -G_fcnp(p,t)).^2)./(G.^2))

where n_exp is the number of points from experimental data.

but, really, Many Thanks again!!!!!

[Star Strider]

As always, my pleasure!

I simply changed the assignments and functions from the provided symbolic code, as best I could.

Even with the revisions, I cannot get it to even come remotely close to fitting ‘G’. If the code is correct, there could be a problems with the constant values or magnitudes.

For the record, the current code is —

G0N=2E5;
Me=18500;
M_star=160000;
k_tdd=1.05E-15;
k_rouse=2E-12;
beta=2.3;
% syms I t mw M0 a b
digits(5)
x = @(t,mw) t./(M_star*k_tdd*(mw.^2));
gx = @(t,mw) -x(t,mw)+sqrt(x(t,mw)).*sqrt(x(t,mw)+(sqrt(pi*x(t,mw)))+pi);
U = @(t,mw) ((1./k_tdd).*(1./(mw.^3)).*t)+(M_star.*gx(t,mw)./mw);
F_sr = @(t,mw) (8/pi^2)*integral(@(I) 1./((2*I-1).^2) .* exp(-(2*I-1).^2 .* U(t,mw)), 1, inf, 'ArrayValued',1);
wGEX = @(mw,M0,a,b) 1.9*(((mw./M0).^(a+1)).*exp(-(mw/M0).^b));
w = @(mw,M0,a,b) wGEX(mw,M0,a,b)./mw;
F_w = @(t,mw,M0,a,b) F_sr(t,mw).*w(mw,M0,a,b);
G = @(t,M0,a,b) integral(@(mw)F_w(t,mw,M0,a,b), 0, inf, 'ArrayValued',1);

G_fcn = @(t,M0,a,b) G0N.*(G(t,M0,a,b)).^beta;

t = [100  50  25  10  5  2.63  1.03  0.526  0.244  0.0991  0.0516  0.0251  0.00996  0.00495  0.00252  0.001 ];

G = [124  879  6110  31900  64000  94400  131000  151000  169000  184000  193000  201000  209000  215000  235000  303000 ];

G_fcnp = @(p,t) G_fcn(t,p(1),p(2),p(3));
 ObjFcn = @(p) norm(G - G_fcnp(p,t))
% n_exp = numel(t);
% ObjFcn = @(p) (1./n_exp).*sum(((G -G_fcnp(p,t)).^2)./(G.^2));
P0=[10 3 200000];
ms = MultiStart;
problem = createOptimProblem('fmincon', 'x0',P0, 'objective',ObjFcn);
P = run(ms,problem,250)
%  [P,resnorm] = lsqcurvefit(G_fcnp, randn(3,1), t, G)

G_out = G_fcnp(P,t)

figure
plot(t, G, 'p')
hold on
plot(t, G_fcnp(P,t), '-r')
hold off
grid

The estimates for the parameters are —

P = 
   336.9743    48.8744   369.3595

It does not converge to a successful result with either version of ‘ObjFcn’.

Answer 2

[3Nz0]

I installed the new version of matlab (2021a), but, again, with this version my results are different from your:

MultiStart completed some of the runs from the start points.

6 out of 250 local solver runs converged with a positive local solver exit flag.

P =

10 3 200000

G_out =

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

I know that it is a very complicated system, for this reason I discretized the problem (file.txt). The problem is I don't know how to extimate the parameters.

[Star Strider]

I have not been able to get a good fit regardless of several different attempts and sevral different approaches, some of which took a few hours to run.

I looked at ‘file.txt’, however I have no idea how to integrate it into the parameter estimation approach.

[3Nz0]

Ok.

Again, thank you so much for your help.

[Star Strider]

As always, my pleasure!

I wish I could have gotten the function to fit the data. Since I do not understand what you are doing or what the function is that you want to fit, unfortunately I cannot offer any other suggestions.

[3Nz0]

Great! I changed the initial values and, finally I have a result (because M0=P(1)).

Thank you for the time you spent due to me.

All the best

[Star Strider]

As always, my pleasure!

What were the succesful parameter estimates?

[3Nz0]

The curve are not completely superimposed but it is fine.

I used fminsearch to evaluate the 3 parameters:

363195.73 26.20 6.13

GREAT!!!

[Star Strider]

Thank you!

I’ll check those. If I can improve on them, I’ll post back here.

[3Nz0]

The results shoul be:

247300; 11.349; 3.366;

I think there is a problem with the constant c=1.9 in wGEX.

This number is a normalizing parameter and the condition is:

c | integral(wGEX)=1.

The final check is not only between G (experimental data) and G_out, but also between wGEX and the experimental data in the attachment.

The .txt is a MWD I digitized from a paper.

I didn't write this just to "simplify" the problem at the beginning