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

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

.

9 个评论
显示 7更早的评论隐藏 7更早的评论

Star Strider 2021-4-30

在 MATLAB Online 中打开

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 2021-4-30

编辑：3Nz0 2021-4-30

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

Star Strider 2021-4-30

编辑：Star Strider 2021-4-30

在 MATLAB Online 中打开

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 2021-5-1

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 2021-5-1

在 MATLAB Online 中打开

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 2021-5-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/818115-how-can-i-use-fmisearch-to-evaluate-the-parameters-in-a-matrix#answer_689370

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.

8 个评论
显示 6更早的评论隐藏 6更早的评论

Star Strider 2021-5-1

Thank you!

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

3Nz0 2021-5-1

MWD.txt

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

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论
显示 7更早的评论隐藏 7更早的评论

更多回答（1 个）

8 个评论
显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论 显示 7更早的评论隐藏 7更早的评论

更多回答（1 个）

8 个评论 显示 6更早的评论隐藏 6更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论