You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to concatenate output of objective function for lsqcurvefit?
1 view (last 30 days)
Show older comments
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 Comments
Hashim
on 21 Nov 2021
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
Accepted Answer
Matt J
on 21 Nov 2021
c0 = [ic_arg{1}(x).'; ic_arg{2}(x).'; ic_arg{3}(x).'];
19 Comments
Hashim
on 21 Nov 2021
Edited: Hashim
on 21 Nov 2021
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
on 21 Nov 2021
You should test your objective function to verify that it works before giving it to lsqcurvefit.
Hashim
on 21 Nov 2021
Edited: Hashim
on 21 Nov 2021
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
on 21 Nov 2021
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
on 21 Nov 2021
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
on 21 Nov 2021
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
on 22 Nov 2021
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
Hashim
on 22 Nov 2021
Then I am sorry I was not able to follow your suggestion regarding defining the s_bulk within the objective function.
Matt J
on 22 Nov 2021
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?
Matt J
on 22 Nov 2021
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
on 23 Nov 2021
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
on 23 Nov 2021
If by rule you mean equations then here is the system of PDEs
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/810404/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/810409/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/810414/image.png)
Results of these PDEs are stored in N*N matrices i.e. (
)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/810419/image.png)
After that I calculate the flux
using pdeval to use in below expression. I use the final row of the solution BTW so basically
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/810424/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/810429/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/810434/image.png)
And thus finally I get my
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/810439/image.png)
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
on 23 Nov 2021
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
on 23 Nov 2021
Edited: Hashim
on 23 Nov 2021
"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
on 23 Nov 2021
Edited: Matt J
on 23 Nov 2021
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
on 23 Nov 2021
Edited: Hashim
on 23 Nov 2021
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?
More Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)