How to concatenate output of objective function for lsqcurvefit?

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);

[Matt J]

We need to see your objective function to know why your efforts are failing.

[Hashim]

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]

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

[Hashim]

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

Answer 1

[Matt J]

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

[Hashim]

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]

You should test your objective function to verify that it works before giving it to lsqcurvefit.

[Hashim]

I can confirm that it works and that I can pass it multiple values via arrayfun and extract multiple outputs like this.

I_anodic       = cell2mat(arrayfun(@PSw_EK_v7, s_bulk, 'un', 0));

I am however unable to do something like this using lsqcurvefit...

lsqcurvefit(cell2mat(arrayfun(@PSw_EK_v7, s_bulk, 'un', 0)), k_m_0, s_bulk, I_exp, lb, ub, options);

I think what I am missing is the inability to output a vector for my array of inputs. Another thing I would like to state here is that for a single input/output value my lsqcurvefit is working i.e. outputting optimized parameters. But it would be much better if lsqcurvefit can take my input array and output the answer vector. I have also tried looping with the unintended side effect that the lsqcurvefit then runs multiple times as well.

[Matt J]

Are you sure it wouldn't be enough just to remove this line

I_anodic    = I_anodic(end);

With arrayfun, you are passing in the same s_bulk multiple times.

[Hashim]

I need the last value from the array not the whole array hence the end. Thing is I need one value of I_anodic against every value of

s_bulk=[0, 5e-06, 10e-06, 20e-06, 50e-06, 100e-06];. 

This is why i've used arrayfun which to my understanding "applies the function func to the elements of A". And so far I am able to generate unique values of I_anodic against s_bulk so I would say that it is working. What is not working is passing the s_bulk to PSw_EK_v7 inside lsqcurvefit to generate unique values of I_anodic. .

[Matt J]

But currently your objective function code uses the entirety of s_bulk to generate one value of I_anodic. If the entirety of s_bulk is available within the workspace of your objective already, you should be able to use it there to calculate all the values you need.

[Hashim]

I tried defining s_bulk inside the function but that did not work. Basically something like this...

s_bulk         = [0, 5e-06, 10e-06, 20e-06, 50e-06, 100e-06];
C.e_layer   = 5.00E-08*ones(length(s_bulk));                          % mol/cm^3
C.m_layer   = 1.3E-06*ones(length(s_bulk));                           % mol/cm^3

[Matt J]

Why did you do that?

[Hashim]

Then I am sorry I was not able to follow your suggestion regarding defining the s_bulk within the objective function.

[Matt J]

No, you are already passing s_bulk to the objective function. My question to you is, why can't you do the entire prediction of I_exp there, in the same call? Why are you returning 1 element from your objective function when you can use s_bulk to return 5?

[Hashim]

Well, that is exactly what I want to know... how do I do that?

[Matt J]

You're the only one who knows your model. Only you understand how the additional elements depend on s_bulk and the unknown parameters.

[Hashim]

The issue is not conceptual rather syntactic. I want to know how can I extract the I_anodic vector against the s_bulk array while using lsqcurvefit?

I have tried.

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

Also,

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

And,

[k_m,resnorm,residual,exitflag,output]=lsqcurvefit(@PSw_EK_v7, k_m_0, [0, 5e-06, 10e-06, 20e-06, 50e-06, 100e-06], I_exp, lb, ub, options);

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.

[Hashim]

If by rule you mean equations then here is the system of PDEs

}

Results of these PDEs are stored in N*N matrices i.e. ()

After that I calculate the flux using pdeval to use in below expression. I use the final row of the solution BTW so basically .

And thus finally I get my

Now my experiemntal data is in 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]

[Matt J]

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]

"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. "

I don't want the predictions written to I_exp (exp is experimental data) rather a new vector on its own I_anodic (pdepe output). I then want to call lsqcurvefit to estimate me the parameter k_m. Now I am able to do that when I pass a single value of s_bulk and k_m so something like this works.

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

And gives me a fitted value of k_m but then the lsqcurvefit is only fitting against 1 output when I want it to fit against the array of I_anodic calculate from each element of s_bulk (one at a time of course).

[Matt J]

My previous comment already tells you how to do that. Here it is again.

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);

I don't want the predictions written to I_exp (exp is experimental data) rather a new vector on its own I_anodic (pdepe output)

Yes, that's why it's refered to as a prediction. I_anodic appears to be your model's prediction of I_exp(i) and depends on s_bulk(i) only, though I've been waiting for you to confirm that.

[Hashim]

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?

[Matt J]

Yes, it will be 6 times slower because lsqcurvefit needs to run the ODE solver for each i=1...6.