第一性原理的最小二乘曲面拟合
Least square surface fitting from first principles
我想做这个 "by hand" 而不是使用曲面拟合工具,因为根据我拥有的数据,曲面拟合可能会有所不同。所以,我首先读取 excel sheet 中的数据,然后初始化一些系数,计算 3D 表面 (f(x,y)) 然后计算总最小二乘和,我会喜欢最小化。每次我 运行 脚本都会告诉我我处于局部最小值,即使我更改了初始值也是如此。更改公差也不影响结果。
这是代码:
% flow function in a separate .m file (approximation, it’s a negative paraboloid, maybe if required, this function may vary):
function Q = flow(P1,P2,a,b,c,d,e,f)
Q1 = a-b.*P1-c.*P1.^2;
Q2 = d-e.*P2-f.*P2.^2;
Q = Q1 + Q2;
% Variable read, I use a xlsread instead
p1a = [-5, -5, -5, -5, -5, -5, -5, -5, -5, -5];
p2a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
qa = [10, 9, 8, 7, 6, 5, 4, 3, 2, 1];
p1b = [-6, -6, -6, -6, -6, -6, -6, -6, -6, -6];
p2b = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
qb = [12, 11, 10, 9, 8, 7, 6, 5, 4, 3];
% Variable initialization
coef = [50, 1, 1, 10, 1, 1];
% Function calculation
q1a = flow(p1a,p2a,coef(1),coef(2),coef(3),coef(4),coef(5),coef(6));
q1b = flow(p1b,p2b,coef(1),coef(2),coef(3),coef(4),coef(5),coef(6));
% Least squares
LQa = (qa-q1a).^2;
LQb = (qb-q1b).^2;
Sa = sum(LQa);
Sb = sum(LQb);
St = Sa+Sb;
% Optimization (minimize the least squares sum)
func = @(coef)(St);
init = coef;
opt = optimoptions('fminunc', 'Algorithm', 'quasi-newton', 'Display', 'iter','TolX', 1e-35, 'TolFun', 1e-30);
[coefmin, Stmin] = fminunc(func, init, opt);
如果你 运行 这个,你应该得到 15546
对于 Stmin
的结果,但是如果你改变系数,你会得到另一个结果,它也会被视为局部最小值。
我做错了什么?
问题是您的 func
只是一个常量。它只是 returns 一个 预先计算的 值 St
,无论您将什么输入传递给 func
,它都是不变的。尝试使用各种不同的输入调用 func
来测试这一点。
您的 objective 函数需要包含让您达到 St
的所有计算。因此,我建议您将 func
替换为保存在 m 文件中的函数,如下所示:
function St = objectiveFunction(coef, p1a, p2a, p1b, p2b, qa, qb, q1a, q1b)
% Function calculation
q1a = flow(p1a,p2a,coef(1),coef(2),coef(3),coef(4),coef(5),coef(6));
q1b = flow(p1b,p2b,coef(1),coef(2),coef(3),coef(4),coef(5),coef(6));
% Least squares
LQa = (qa-q1a).^2;
LQb = (qb-q1b).^2;
Sa = sum(LQa);
Sb = sum(LQb);
St = Sa+Sb;
end
然后在您的脚本调用 objectiveFunction
中使用这样的匿名函数:
[coefmin, Stmin] = fminunc(@(coef)(objectiveFunction(coef, p1a, p2a, p1b, p2b, qa, qb, q1a, q1b)), init, opt);
这个想法是创建一个匿名函数,它只接受一个参数,coef
,这是 fminunc
将扰乱并传递回您的 objective 函数的变量。您的 objectiveFunction
需要的其他参数(即 p1a, p2a, p1b,...
)现在被认为是由您的匿名函数预先计算的,因此由 fminunc
.
其余代码可以保持不变。
我想做这个 "by hand" 而不是使用曲面拟合工具,因为根据我拥有的数据,曲面拟合可能会有所不同。所以,我首先读取 excel sheet 中的数据,然后初始化一些系数,计算 3D 表面 (f(x,y)) 然后计算总最小二乘和,我会喜欢最小化。每次我 运行 脚本都会告诉我我处于局部最小值,即使我更改了初始值也是如此。更改公差也不影响结果。
这是代码:
% flow function in a separate .m file (approximation, it’s a negative paraboloid, maybe if required, this function may vary):
function Q = flow(P1,P2,a,b,c,d,e,f)
Q1 = a-b.*P1-c.*P1.^2;
Q2 = d-e.*P2-f.*P2.^2;
Q = Q1 + Q2;
% Variable read, I use a xlsread instead
p1a = [-5, -5, -5, -5, -5, -5, -5, -5, -5, -5];
p2a = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
qa = [10, 9, 8, 7, 6, 5, 4, 3, 2, 1];
p1b = [-6, -6, -6, -6, -6, -6, -6, -6, -6, -6];
p2b = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
qb = [12, 11, 10, 9, 8, 7, 6, 5, 4, 3];
% Variable initialization
coef = [50, 1, 1, 10, 1, 1];
% Function calculation
q1a = flow(p1a,p2a,coef(1),coef(2),coef(3),coef(4),coef(5),coef(6));
q1b = flow(p1b,p2b,coef(1),coef(2),coef(3),coef(4),coef(5),coef(6));
% Least squares
LQa = (qa-q1a).^2;
LQb = (qb-q1b).^2;
Sa = sum(LQa);
Sb = sum(LQb);
St = Sa+Sb;
% Optimization (minimize the least squares sum)
func = @(coef)(St);
init = coef;
opt = optimoptions('fminunc', 'Algorithm', 'quasi-newton', 'Display', 'iter','TolX', 1e-35, 'TolFun', 1e-30);
[coefmin, Stmin] = fminunc(func, init, opt);
如果你 运行 这个,你应该得到 15546
对于 Stmin
的结果,但是如果你改变系数,你会得到另一个结果,它也会被视为局部最小值。
我做错了什么?
问题是您的 func
只是一个常量。它只是 returns 一个 预先计算的 值 St
,无论您将什么输入传递给 func
,它都是不变的。尝试使用各种不同的输入调用 func
来测试这一点。
您的 objective 函数需要包含让您达到 St
的所有计算。因此,我建议您将 func
替换为保存在 m 文件中的函数,如下所示:
function St = objectiveFunction(coef, p1a, p2a, p1b, p2b, qa, qb, q1a, q1b)
% Function calculation
q1a = flow(p1a,p2a,coef(1),coef(2),coef(3),coef(4),coef(5),coef(6));
q1b = flow(p1b,p2b,coef(1),coef(2),coef(3),coef(4),coef(5),coef(6));
% Least squares
LQa = (qa-q1a).^2;
LQb = (qb-q1b).^2;
Sa = sum(LQa);
Sb = sum(LQb);
St = Sa+Sb;
end
然后在您的脚本调用 objectiveFunction
中使用这样的匿名函数:
[coefmin, Stmin] = fminunc(@(coef)(objectiveFunction(coef, p1a, p2a, p1b, p2b, qa, qb, q1a, q1b)), init, opt);
这个想法是创建一个匿名函数,它只接受一个参数,coef
,这是 fminunc
将扰乱并传递回您的 objective 函数的变量。您的 objectiveFunction
需要的其他参数(即 p1a, p2a, p1b,...
)现在被认为是由您的匿名函数预先计算的,因此由 fminunc
.
其余代码可以保持不变。