Can lsqnonlin regress morethan one parameter?
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Hi everyone,
I have a code that regresses one parameter well (see the code below):
%% Regression Example
%% Parameters
p.KW = 10^-7;
p.K1 = 10^-6.35;
p.K2 = 10^-10.3;
%% Initial conditions
y0.cB0 = 3.5;
y0.cT0 = 8;
%% Known values
p.cNa = 2*y0.cT0;
p.cCl = 12;
%% Measurements
t = [0 2 4 6 8 10];
Data.t = [0 2 4 6 8 10];
Data.cB_measured_1 = [3.76883357 2.143052955 -0.589243432 0.689078443 0.286944142 -0.589693715];
Data.cB_measured_2 = [3.283203989 1.397422681 2.329378469 1.642710298 -0.547381948 1.581612167];
Data.cB_measured_3 = [3.862702112 1.194583011 0.897551451 0.155508754 0.065489348 0.808999237];
Data.pH_measured_1 = [7.006085712 6.905004758 6.85705829 6.819587483 6.829728045 6.834313912];
Data.pH_measured_2 = [6.996884305 6.901179764 6.85761217 6.828627943 6.825494374 6.810138731];
Data.pH_measured_3 = [7.000879323 6.879362133 6.839654615 6.823567365 6.793112817 6.832395362];
%% Regressing using lsqnonlin,
%[p.RegrPara1, p.RegrPara2] = lsqnonlin(@(k) funcLSQ(k, Data, p, y0), 1, 1e-4, 10);
p.RegrPara1 = lsqnonlin(@(k) funcLSQ(k, Data, p, y0), 1, 1e-4, 10);
[cB, cT, pH] = funcSimulate(Data.t, y0.cB0, y0.cT0, p);
disp(['RegrPara1 = ',num2str(p.RegrPara1)]);
%% Plotting the results
subplot(2,1,1)
plot(t, cB, t, cT, t, Data.cB_measured_1,'ko',t, Data.cB_measured_2,'ko',t, Data.cB_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('Concentration'); legend('C_B','C_T','Measured C_B');
subplot(2,1,2)
plot(t, pH, t, Data.pH_measured_1,'ko',t, Data.pH_measured_2,'ko',t, Data.pH_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('pH'); legend('pH','Measured pH');
%-------------------------------------------------------------------------------------------------------------------------------
%% Function used to calculate the error vector (which will be used by lsqnonlin)
function Err = funcLSQ(RegrPara1, Data, p, y0)
p.RegrPara1 = RegrPara1;
% p.RegrPara2 = RegrPara2;
t = [0 2 4 6 8 10];
% Using the guessed value of "k" calculate cB and pH
[cB, ~, pH] = funcSimulate(t, y0.cB0, y0.cT0, p);
% Difference between measured and predicted cB
Err_cB = [ Data.cB_measured_1 - cB; Data.cB_measured_2 - cB; Data.cB_measured_3 - cB];
% Difference between measured and predicted pH
Err_pH = [Data.pH_measured_1 - pH; Data.pH_measured_2 - pH; Data.pH_measured_3 - pH];
% Error vector containing both the cB and pH error.
Err = [Err_cB; Err_pH];
end
%% Calcuating state and non-state variables
function [cB, cT, pH] = funcSimulate(t, cB0, cT0, p)
[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cB0; cT0]);
cB = y(:,1);
cT = y(:,2);
pH = zeros(length(t), 1);
for i = 1 : length(t)
[~, pH(i)] = funcOtherVariables(cB(i), cT(i), p);
end
end
%% ODES
function dy_dt = funcODE(t, y, p)
cB = y(1);
cT = y(2);
r = funcOtherVariables(cB, cT, p);
dcB_dt = -r;
dcT_dt = -r;
dy_dt = [dcB_dt; dcT_dt];
end
%% Other variables
function [r, pH] = funcOtherVariables(cB, cT, p)
% Determines pH by solving the charge balance
pH = fzero(@(pH) funcCharge(pH, cB, cT, p), 7);
% Determines the reaction rate once the pH (and thus cHCO3) is known
cH = 10^-pH;
cHCO3 = ( p.K1/cH / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT;
r(1) = p.RegrPara1*cB*cHCO3;
% r(2) = p.RegrPara2*cB*cHCO3;
% r = [r(1) r(2)];
r = r(1);
end
%% Solving algebraic equation using fsolve
function z = funcCharge(pH, cB, cT, p)
cH = 10^-pH;
% Left-hand side of the charge balance
LHS = cB + p.cNa + cH;
% Right-hand side of the charge balance
RHS = p.KW/cH + ( (p.K1/cH + 2*p.K1*p.K2/cH^2) / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT + p.cCl;
z = LHS - RHS; % Should be equal to zero
end
However, when I introduce more paraters, it gives error message (also see code below):
%% Regression Example
%% Parameters
p.KW = 10^-7;
p.K1 = 10^-6.35;
p.K2 = 10^-10.3;
%% Initial conditions
y0.cB0 = 3.5;
y0.cT0 = 8;
%% Known values
p.cNa = 2*y0.cT0;
p.cCl = 12;
%% Measurements
t = [0 2 4 6 8 10];
Data.t = [0 2 4 6 8 10];
Data.cB_measured_1 = [3.76883357 2.143052955 -0.589243432 0.689078443 0.286944142 -0.589693715];
Data.cB_measured_2 = [3.283203989 1.397422681 2.329378469 1.642710298 -0.547381948 1.581612167];
Data.cB_measured_3 = [3.862702112 1.194583011 0.897551451 0.155508754 0.065489348 0.808999237];
Data.pH_measured_1 = [7.006085712 6.905004758 6.85705829 6.819587483 6.829728045 6.834313912];
Data.pH_measured_2 = [6.996884305 6.901179764 6.85761217 6.828627943 6.825494374 6.810138731];
Data.pH_measured_3 = [7.000879323 6.879362133 6.839654615 6.823567365 6.793112817 6.832395362];
%% Regressing using lsqnonlin,
[p.RegrPara1, p.RegrPara2] = lsqnonlin(@(RegrPara1, RegrPara2) funcLSQ(RegrPara1, RegrPara2, Data, p, y0), 1, 1e-4, 10);
[cB, cT, pH] = funcSimulate(Data.t, y0.cB0, y0.cT0, p);
disp(['k = ',num2str(p.k)]);
%% Plotting the results
subplot(2,1,1)
plot(t, cB, t, cT, t, Data.cB_measured_1,'ko',t, Data.cB_measured_2,'ko',t, Data.cB_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('Concentration'); legend('C_B','C_T','Measured C_B');
subplot(2,1,2)
plot(t, pH, t, Data.pH_measured_1,'ko',t, Data.pH_measured_2,'ko',t, Data.pH_measured_3,'ko','LineWidth',1);
xlabel('Time'); ylabel('pH'); legend('pH','Measured pH');
%-------------------------------------------------------------------------------------------------------------------------------
%% Function used to calculate the error vector (which will be used by lsqnonlin)
function Err = funcLSQ(t, RegrPara1,RegrPara2, Data, p, y0)
p.RegrPara1 = RegrPara1;
p.RegrPara2 = RegrPara1;
t = [0 2 4 6 8 10];
% Using the guessed value of "k" calculate cB and pH
[cB, cT, pH] = funcSimulate(t, y0.cB0, y0.cT0, p);
% Difference between measured and predicted cB
Err_cB = [ Data.cB_measured_1 - cB; Data.cB_measured_2 - cB; Data.cB_measured_3 - cB];
% Difference between measured and predicted pH
Err_pH = [Data.pH_measured_1 - pH; Data.pH_measured_2 - pH; Data.pH_measured_3 - pH];
% Error vector containing both the cB and pH error.
Err = [Err_cB; Err_pH];
end
%% Calcuating state and non-state variables
function [cB, cT, pH] = funcSimulate(t, cB0, cT0, p)
[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cB0; cT0]);
cB = y(:,1);
cT = y(:,2);
pH = zeros(length(t), 1);
for i = 1 : length(t)
[~, pH(i)] = funcOtherVariables(cB(i), cT(i), p);
end
end
%% ODES
function dy_dt = funcODE(t, y, p)
cB = y(1);
cT = y(2);
r = funcOtherVariables(cB, cT, p);
dcB_dt = -r;
dcT_dt = -r;
dy_dt = [dcB_dt; dcT_dt];
end
%% Other variables
function [r, pH] = funcOtherVariables(cB, cT, p)
% Determines pH by solving the charge balance
pH = fzero(@(pH) funcCharge(pH, cB, cT, p), 7);
% Determines the reaction rate once the pH (and thus cHCO3) is known
cH = 10^-pH;
cHCO3 = ( p.K1/cH / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT;
r(1) = p.RegrPara1*cB*cHCO3;
r(2) = p.RegrPara2*cB*cHCO3;
r = [r(1) r(2)];
end
%% Solving algebraic equation using fsolve
function z = funcCharge(pH, cB, cT, p)
cH = 10^-pH;
% Left-hand side of the charge balance
LHS = cB + p.cNa + cH;
% Right-hand side of the charge balance
RHS = p.KW/cH + ( (p.K1/cH + 2*p.K1*p.K2/cH^2) / (1 + p.K1/cH + p.K1*p.K2/cH^2) )*cT + p.cCl;
z = LHS - RHS; % Should be equal to zero
end
The error message is :
Not enough input arguments.
Error in RegressionExample>@(RegrPara1,RegrPara2)funcLSQ(RegrPara1,RegrPara2,Data,p,y0) (line 29)
[p.RegrPara1, p.RegrPara2] = lsqnonlin(@(RegrPara1, RegrPara2) funcLSQ(RegrPara1, RegrPara2, Data, p, y0), 1, 1e-4, 10);
Error in lsqnonlin (line 218)
initVals.F = feval(funfcn{3},xCurrent,varargin{:});
Error in RegressionExample (line 29)
[p.RegrPara1, p.RegrPara2] = lsqnonlin(@(RegrPara1, RegrPara2) funcLSQ(RegrPara1, RegrPara2, Data, p, y0), 1, 1e-4, 10);
Caused by:
Failure in initial objective function evaluation. LSQNONLIN cannot continue.
Hence my question is, can lsqnonlin regress morethan one parameter or not?
0 comentarios
Respuestas (1)
Matt J
el 30 de Jun. de 2021
Editada: Matt J
el 30 de Jun. de 2021
It is a strange usage of lsqnonlin to solve for only a single unknown parameter. For that, you would normally just use fminsearch or fminbnd. In the more usual case where there is more than a single unknown, lsqnonlin expects your objective function to accept it in the form of a vector:
objective = @(x) funcLSQ(x(1), x(2), Data, p, y0)
11 comentarios
Matt J
el 30 de Jun. de 2021
The differences between lsqnonlin and lsqcurvefit are very minor. Basically, lsqcurvefit is more convenient if you want to plot a fitted curve after the optimization.
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!