Can lsqnonlin regress morethan one parameter?
Show older comments
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?
Answers (1)
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 Comments
Dursman Mchabe
on 30 Jun 2021
Matt J
on 30 Jun 2021
You have given additional arguments to lsqnonlin (1,1e-4,10) which for some reason are all scalars. They should not be.The second argument to lsqnonlin (your initial guess) also needs to be a vector containing all the unknowns. The same for any upper and lower bound (lb and ub) arguments you wish to give.
Dursman Mchabe
on 30 Jun 2021
Matt J
on 30 Jun 2021
A few problems.
First, you have given as your initial guess a 3-vector RegrPara0=[1, 1e-4, 10]. This tells lsqnonlin, among other things, that your objective function will accept 3 unknowns. However, in the expression
@(RegrPara) funcLSQ(RegrPara(1), RegrPara(2), Data, p,[1 1])
it appears that RegrPara contains only 2 unknown elements RegrPara(1) and RegrPara(2).
Additionally, you are passing 5 arguments to funcLSQ, but you have written it to expect 6:
function Err = funcLSQ(t, RegrPara1,RegrPara2, Data, p, y0)
Dursman Mchabe
on 30 Jun 2021
Matt J
on 30 Jun 2021
Now you have written funcLSQ to accept 4 arguments but you are passing in 5...
Dursman Mchabe
on 30 Jun 2021
Matt J
on 30 Jun 2021
Well, before you do anything with lsqnonlin, I recommend that you first make sure that funcLSQ is working. When I simply call it with your initial point RegPara=[1e-4,10],
funcLSQ([1e-4,10], Data, p,y0);
I obtain the following errors
Error using odearguments (line 93)
@(T,Y)FUNCODE(T,Y,P) must return a column vector.
Error in ode23s (line 121)
= odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in test>funcSimulate (line 59)
[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cB0; cT0]);
Error in test>funcLSQ (line 43)
[cB, cT, pH] = funcSimulate(t, y0.cB0, y0.cT0, p);
Error in test (line 30)
funcLSQ([1e-4,10], Data, p,y0);
Dursman Mchabe
on 30 Jun 2021
Dursman Mchabe
on 30 Jun 2021
Matt J
on 30 Jun 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.
Categories
Find more on Support Vector Machine Regression in Help Center and File Exchange
on 29 Jun 2021
on 30 Jun 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
