現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
How to concatenate output of objective function for lsqcurvefit?
1 回表示 (過去 30 日間)
古いコメントを表示
Hashim
2021 年 11 月 21 日
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 件のコメント
Hashim
2021 年 11 月 21 日
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 日
c0 = [ic_arg{1}(x).'; ic_arg{2}(x).'; ic_arg{3}(x).'];
19 件のコメント
Hashim
2021 年 11 月 21 日
編集済み: Hashim
2021 年 11 月 21 日
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 月 21 日
You should test your objective function to verify that it works before giving it to lsqcurvefit.
Hashim
2021 年 11 月 21 日
編集済み: Hashim
2021 年 11 月 21 日
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
2021 年 11 月 21 日
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
2021 年 11 月 21 日
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
2021 年 11 月 21 日
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
2021 年 11 月 22 日
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
2021 年 11 月 22 日
Then I am sorry I was not able to follow your suggestion regarding defining the s_bulk within the objective function.
Matt J
2021 年 11 月 22 日
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
2021 年 11 月 22 日
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
2021 年 11 月 23 日
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
2021 年 11 月 23 日
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
2021 年 11 月 23 日
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 日
"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
2021 年 11 月 23 日
編集済み: Matt J
2021 年 11 月 23 日
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
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?
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- 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)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)