How to concatenate output of objective function for lsqcurvefit?

Question

Hashim 2021-11-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1591959-how-to-concatenate-output-of-objective-function-for-lsqcurvefit

编辑： Hashim 2021-11-23

采纳的回答： Matt J

在 MATLAB Online 中打开

Hi All,

I have an objective function of the form

function [I_num] = PSw_EK_v7(k_m, s_bulk)

k_m is the variable I want to optimize while s_bulk is an array of input data which I want to feed the input function. I am using the lsqcurvefit/lsqnonlin function to fit my simulation to experimental data of the form.

s_bulk         = [0, 5e-06, 10e-06, 20e-06, 50e-06, 100e-06];
I_exp          = [0, 2.35e-03, 2.7e-03, 3.125e-03, 3.29e-03, 3.25e-03]

Now what I am trying to do is feed the lsqcurvefit function my s_bulk array and get an array of I_num values but I am having difficulty doing that. I have tried arrayfun function but that does not seem to work inside the lsqcurvefit function.

So, to put things short I want to give my s_bulk array and get an array of I_num because that is what lsqcurvefit expects.

[k_m,resnorm,residual,exitflag,output]=lsqcurvefit(@PSw_EK_v7, k_m_0, s_bulk, I_exp, lb, ub, options);

4 个评论
显示 2更早的评论隐藏 2更早的评论

Hashim 2021-11-21

在 MATLAB Online 中打开

function [I_anodic] = PSw_EK_v7(k_m, s_bulk)
% Constants Used in the Nernst Equation   
C.n         = 1;                                 % No. of Electrons Involved
C.F         = 96485.3329;                        % C/mol
C.R         = 8.314;                             % J/(mol.K) or CV/(mol.K)
C.Temp      = 273.15+37;                         % K
% Diffusivity Constants 
C.D_s       = 7.00E-06;                          % cm^2/s
C.D_m       = 2.81E-06;                          % cm^2/s (De of polymer)
C.D_e       = 0;                                 % because covalently bound
% k_e = kcat / Km + S_bulk hence parameter above is unused
% C.k_m       = 4.92E+07;                          % (mol/cm3)^{-1}s^{-1}
C.kcat      = 800;                               % Turnover No. (1/s) 
C.k_e       = C.kcat;                            % Refer to Schnell
C.KM        = 2.00E-05;                          % Michalis Constant (mol/cm^3)
C.P_s       = 1;                                 % Partition Coefficient
C.P_m       = 1;                                 % Partition Coefficient
% Initial Concentrations 
C.e_layer   = 5.00E-08;                          % mol/cm^3
C.m_layer   = 1.3E-06;                           % mol/cm^3
C.Area      = 0.0707;                            % cm^2 
C.l         = 0.001;                             % cm 
m           = 0;                                 % Cartesian Co-ordinates
N           = 20;                                % No. of Points
 
C.tspan1    = linspace(0, 500, N);               % s 
xmesh       = linspace(0, C.l, N);               % cm 
% xmesh       = logspace(log10(0.000006), log10(C.l), N); xmesh(1) = 0;
% FIRST LINEAR POTENTIAL SWEEP PDEPE Solver Command
E           = 0.0;                               % V
C.E0        = 0.2;                               % V
C.ScanRt    = 0.001;                             % V/s
C.E_1       = E+(C.ScanRt.*C.tspan1);            % Potential Sweep 1
C.epsilon1  = ((C.n.*C.F)./(C.R.*C.Temp)).*(C.E_1-C.E0); 
C.c_mred1   = C.m_layer./(1+exp(C.epsilon1));    % Mred
% FIRST POTENTIAL SWEEP Initial & Boundary Condition Call
ic_arg      = {@(x)s_bulk.*ones(size(N)) ; @(x)C.e_layer.*ones(size(N)); ...
               @(x)C.m_layer.*ones(size(N))}; 
IC          =  @(x)PDE_PSw_EK_IC(x, ic_arg, k_m); 
BC          =  @(xl, cl, xr, cr, t)PDE_PSw_EK_BC(xl, cl, xr, cr, t, C, s_bulk, k_m); 
optns      = odeset('MaxStep',1e-01,'RelTol',1e-09,'AbsTol',1e-10); 
sol1       = pdepe(m, @(x, t, c, DcDx)PDE_PSw_EK(x, t, c, DcDx, C, k_m), ... 
                   IC, BC, xmesh, C.tspan1, optns);
               
% Concentration Profiles c(i, j, k)(Solutions) 
c1          = sol1(:, :, 1);                     % Substrate Conc.  
c2          = sol1(:, :, 2);                     % Ered Conc.
c3          = sol1(:, :, 3);                     % Mred Conc.
c4          = C.m_layer-c3;                      % Mox Conc. 
% Calculation of Mox flux at electrode surface 
for counter = 2:N
    [~, dc4Mdx_0(counter)]  = pdeval(m, xmesh, c4(counter,:),  0);  % Mox Flux at x=0
end 
j_num_ano   = (-C.D_m .* dc4Mdx_0);
I_anodic    = (C.n .* C.F .* j_num_ano)./(C.Area); 
I_anodic    = I_anodic(end);
z = 1; 
%% pdepe Function PSw_EK_v8
%%
function [cc, ff, ss] = PDE_PSw_EK(x, t, c, DcDx, C, k_m)
%          Substrate;  Ered;  Mred;   
cc          = [    1;     1;     1]; 
ff          = [C.D_s; C.D_e; C.D_m].*DcDx;
menten      = (C.k_e.*c(1)).*(C.e_layer-c(2))./(C.KM+c(1));
mred        = k_m.*(C.m_layer-c(3)).*c(2);
ss          = [-menten; -mred+menten; +mred]; 
           
end
%% Initial Condition --> Initial Concentrations of Species 
%%
function c0 = PDE_PSw_EK_IC(x, ic_arg, k_m)
   
%                 Substrate;         Ered;          Mred;   
c0          = [ic_arg{1}(x); ic_arg{2}(x); ic_arg{3}(x)];
end
%% Boundary Condition - 1 
%%
function [pl, ql, pr, qr] = PDE_PSw_EK_BC(xl, cl, xr, cr, t, C, s_bulk, k_m)
% ---------------------
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% -------|-------------
%        pl           pr 
% Order --> Substrate; Mox; Mred; Eox; Ered 
%              Substrate;Er;                                  Mox;   
pl     = [0             ; 0; cl(3)-interp1(C.tspan1,C.c_mred1,t)];  
 
ql     = [1             ; 1; 0];
pr     = [cr(1)-s_bulk  ; 0; 0];
qr     = [0             ; 1; 1]; 
end
end

Matt J 2021-11-21

编辑：Matt J 2021-11-21

Also, can you provide to us k_m_0, lb, and ub?

Hashim 2021-11-21

在 MATLAB Online 中打开

lb             = 1e06;
ub             = 1e09; 
k_m_0          = 1.9e07;

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2021-11-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1591959-how-to-concatenate-output-of-objective-function-for-lsqcurvefit#answer_836654

在 MATLAB Online 中打开

c0 = [ic_arg{1}(x).'; ic_arg{2}(x).'; ic_arg{3}(x).'];

19 个评论
显示 17更早的评论隐藏 17更早的评论

Hashim 2021-11-21

编辑：Hashim 2021-11-21

在 MATLAB Online 中打开

Ran it with the following commands.

s_bulk         = [0, 5e-06];
I_exp          = [0, 2.35e-030];
[k_m,resnorm,residual,exitflag, output]=lsqcurvefit(@PSw_EK_v7, k_m_0, s_bulk, I_exp, lb, ub, options);

Error message

Matrix dimensions must agree.
Error in PSw_EK_v7/PDE_PSw_EK (line 100)
ff          = [C.D_s; C.D_e; C.D_m].*DcDx;
Error in PSw_EK_v7>@(x,t,c,DcDx)PDE_PSw_EK(x,t,c,DcDx,C,k_m) (line 64)
sol1       = pdepe(m, @(x, t, c, DcDx)PDE_PSw_EK(x, t, c, DcDx, C, k_m), ...
Error in pdepe (line 246)
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});
Error in PSw_EK_v7 (line 64)
sol1       = pdepe(m, @(x, t, c, DcDx)PDE_PSw_EK(x, t, c, DcDx, C, k_m), ...
Error in lsqcurvefit (line 202)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});
Error in PSw_EK_v7_Caller (line 24)
    lsqcurvefit(@PSw_EK_v7, k_m_0, s_bulk, I_exp, lb, ub, options);
Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

Matt J 2021-11-23

在 MATLAB Online 中打开

I just want to pass the s_bulk array to lsqcurvefit so that I don't have to run for every individual element of s_bulk.

You are already doing that. In your original post, you are calling with the syntax,

[k_m,resnorm,residual,exitflag,output]=...
    lsqcurvefit(@PSw_EK_v7, k_m_0, s_bulk, I_exp, lb, ub, options);

With this syntax, the entire vector s_bulk is being passed to PSw_EK_v7 every time lsqcurvefit calls it. You can verify this in many ways, but the simplest would be to add a line at the top of PSw_EK_v7 that prints s_bulk to the screen

function [I_anodic] = PSw_EK_v7(k_m, s_bulk)
  
 display(s_bulk)
 
...
end

If you do this, you will see that the entire 1x6 vector s_bulk is available to PSw_EK_v7() every time it is called.

Because s_bulk is being passed to PSw_EK_v7 as a vector, it is not clear why you haven't written PSw_EK_v7 to return a 1x6 vector of predictions to I_exp. You have everything you need within the workspace of PSw_EK_v7.

It is also not clear what the dependence should be. You seem to be suggesting that each element I_exp(i) depends only on a single corresponding s_bulk(i) and that PSw_EK_v7 has been written to expect only a single s_bulk(i) as input. If so, that is clear to no one but you. As far as we and lsqcurvefit know, your model function PSw_EK_v7 is an arbitrary mapping from R^6 to R^6. There is not even a requirement that s_bulk and I_exp be the same size, in general. s_bulk could have been a 1000x1000 matrix even if I_exp is still a 1x6 vector and lsqcurvefit would not have cared.

But if it is true that PSw_EK_v7 has been written to expect only a single s_bulk(i) as input, you could do this instead,

wrapper=@(k_m,s_bulk) arrayfun(@(si)PSw_EK_v7(k_m,si), s_bulk);
[k_m,resnorm,residual,exitflag,output]=...
lsqcurvefit(wrapper, k_m_0, s_bulk, I_exp, lb, ub, options);

However, it may be more efficient to vectorize the input to the ODE solver in some way...

Hashim 2021-11-23

编辑：Hashim 2021-11-23

Hi, it might be working, let me confirm it and get it back to you. Apparently the execution time is increased manyfolds when I run it this way. And please accept my humble gratitude for your help so far. Is there any way I can get the I_anodic as an output while calling lsqcurvefit?