3D 中的 MATLAB 曲线拟合，具有附加边界
MATLAB Curve Fitting in 3D, with additional boundaries
简介
假设我有一组实验数据，需要找到描述所选系列的多项式近似值。实验结果取决于两个变量——时间和浓度。让示例数据如下所示：
Experiment=[1.5 0.2 0.4 0.4 0.2 0.2 2.0 0.2 0.4 0.4 0.2 0.2];
Time=[0 5 10 0 5 10 0 5 10 0 5 10];
Concentration=[0 0 0 1 1 1 2 2 2 3 3 3];
多项式可以很容易地拟合和绘制如下：
Time = transpose(Time);
Concentration = transpose(Concentration);
Experiment= transpose(Experiment);
f1 = fit( [Time, Concentration], Experiment, 'poly23' );
pl=plot(f1, [Time, Concentration], Experiment);
问题
上面描述的简单程序完全没问题，并给出了一个多项式图：
但当时间为4-10且浓度小于1时，多项式结果为负。我正在研究的系统是生物系统。所以任何负值在物理上都是不可能的。所以我的问题是：如何设置任何 boundaries/constraints 以防止结果多项式在实验范围内为负数？ 如何强制 MATLAB 给我近似值，这给出 Z>如果时间在 0 到 10 之间，浓度在 0 到 3 之间，则始终为 0?
对于使用 fmincon 的非线性约束优化，您首先需要定义一个确定 z 的函数（即预测 x 和 y 的结果）：
function z = poly_model(x, y, p)
% extract parameters
p00 = p(1);
p10 = p(2);
p01 = p(3);
p20 = p(4);
p11 = p(5);
p02 = p(6);
p21 = p(7);
p12 = p(8);
p03 = p(9);
% poly23 model
z = p00 + p10 .* x + p01 .* y + p20 .* x.^2 + p11 .* x .* y + ...
p02 .* y.^2 + p21 .* x.^2 .* y + p12 .* x .* y.^2 + p03 .* y.^3;
end
请注意，所有乘法和幂都是按元素计算的（.* 和 .^）。这允许评估 x 和 y 的矩阵输入函数，这是计算要在实验数据范围内施加的约束所必需的。
约束已在单独的函数中定义。来自文档：
Nonlinear constraints, specified as a function handle or function name. nonlcon is a function that accepts a vector or array x and returns two arrays, c(x) and ceq(x).
c(x) is the array of nonlinear inequality constraints at x. fmincon attempts to satisfy
c(x) <= 0 for all entries of c.
ceq(x) is the array of nonlinear equality constraints at x. fmincon attempts to satisfy
ceq(x) = 0 for all entries of ceq.
So in your case, the constraint function can be defined as:
function [c, ceq] = constraint_eq(x, y, p)
% evaluate the model for required x and y
z_model = poly_model(x, y, p);
% and constrain z to be positive:
c = -z_model; % z_model >= 0, c(p) <= 0, hence c = -z_model
% no inequality constraint needed
ceq = [];
end
接下来，您需要定义一个优化函数，使实验数据与模型预测之间的误差最小化：
function err = cost_function(x, y, z, p)
z_model = poly_model(x, y, p); % determine model prediction z for x and y
ev = z_model - z; % error vector
err = norm(ev, 2)^2; % sum of squared error
end
最后，调用 fmincon 例程：
clc
clear
close all
% data
Experiment = [1.5 0.2 0.4 0.4 0.2 0.2 2.0 0.2 0.4 0.4 0.2 0.2];
Time = [0 5 10 0 5 10 0 5 10 0 5 10];
Concentration = [0 0 0 1 1 1 2 2 2 3 3 3];
% short notation for readability
x = Time;
y = Concentration;
z = Experiment;
% define XV and YV to fulfil constraint over the entire x and y domain
xv = linspace(min(x), max(x), 20);
yv = linspace(min(y), max(y), 20);
[XV, YV] = meshgrid(xv, yv);
% initial guess parameters?
p0 = rand(9, 1);
p_final = fmincon(@(p) cost_function(x, y, z, p), p0, [], [], [], [], [], [], @(p) constraint_eq(XV, YV, p));
%% check result:
ZV = poly_model(XV, YV, p_final); % evaluate points in poly23 plane
% plot result
figure(1); clf;
scatter3(x, y, z, 200, 'b.');
hold on;
surf(XV, YV, ZV)
初始参数的影响p0
正如@James Philips 在评论中指出的那样，您还可以使用无约束优化的解决方案作为约束优化的起点。对于提供的实验数据和选择的模型，你会发现实际上并没有什么区别：
% The random initial guess:
p0 = rand(9, 1);
% Optimal solution for random p0
p_rand = fmincon(@(p) cost_function(x, y, z, p), p0, [], [], [], [], [], [], @(p) constraint_eq(XV, YV, p));
% first running unconstrained optimization and use p_unc
% as start point for constrained optimization
p_unc = fmincon(@(p) cost_function(x, y, z, p), p0, [], []);
p_con= fmincon(@(p) cost_function(x, y, z, p), p_unc, [], [], [], [], [], [], @(p) constraint_eq(XV, YV, p));
% Compare errors:
SSE_unc = cost_function(x,y,z,p_unc)
SSE_con = cost_function(x,y,z,p_con)
SSE_rand = cost_function(x,y,z,p_rand)
% compare poly23 parameters
p_all = [p_unc, p_con, p_rand]
这将给出：
SSE_unc =
1.0348
SSE_con =
1.1889
SSE_rand =
1.1889
p_all =
1.3375 1.2649 1.2652
-0.3425 -0.2617 -0.2618
-1.6069 -1.0620 -1.0625
0.0258 0.0187 0.0187
0.0175 -0.0018 -0.0016
1.5708 1.0717 1.0721
-0.0042 -0.0018 -0.0018
0.0125 0.0094 0.0094
-0.3722 -0.2627 -0.2628
在这种情况下，您会发现所找到的参数存在非常小的差异，但求解器很可能需要较少的迭代次数才能获得此解。通过调整求解器设置（最优容差和约束容差），p_rand 和 p_con 的解会更接近。
一般来说，检查多个随机初始猜测是一种很好的做法，以确保您没有找到局部最小值（例如通过使用 MultiStart）。
简介
假设我有一组实验数据，需要找到描述所选系列的多项式近似值。实验结果取决于两个变量——时间和浓度。让示例数据如下所示：
Experiment=[1.5 0.2 0.4 0.4 0.2 0.2 2.0 0.2 0.4 0.4 0.2 0.2];
Time=[0 5 10 0 5 10 0 5 10 0 5 10];
Concentration=[0 0 0 1 1 1 2 2 2 3 3 3];
多项式可以很容易地拟合和绘制如下：
Time = transpose(Time);
Concentration = transpose(Concentration);
Experiment= transpose(Experiment);
f1 = fit( [Time, Concentration], Experiment, 'poly23' );
pl=plot(f1, [Time, Concentration], Experiment);
问题
上面描述的简单程序完全没问题，并给出了一个多项式图：
对于使用 fmincon 的非线性约束优化，您首先需要定义一个确定 z 的函数（即预测 x 和 y 的结果）：
function z = poly_model(x, y, p)
% extract parameters
p00 = p(1);
p10 = p(2);
p01 = p(3);
p20 = p(4);
p11 = p(5);
p02 = p(6);
p21 = p(7);
p12 = p(8);
p03 = p(9);
% poly23 model
z = p00 + p10 .* x + p01 .* y + p20 .* x.^2 + p11 .* x .* y + ...
p02 .* y.^2 + p21 .* x.^2 .* y + p12 .* x .* y.^2 + p03 .* y.^3;
end
请注意，所有乘法和幂都是按元素计算的（.* 和 .^）。这允许评估 x 和 y 的矩阵输入函数，这是计算要在实验数据范围内施加的约束所必需的。
约束已在单独的函数中定义。来自文档：
Nonlinear constraints, specified as a function handle or function name. nonlcon is a function that accepts a vector or array x and returns two arrays, c(x) and ceq(x).
c(x) is the array of nonlinear inequality constraints at x. fmincon attempts to satisfy
c(x) <= 0 for all entries of c.
ceq(x) is the array of nonlinear equality constraints at x. fmincon attempts to satisfy
ceq(x) = 0 for all entries of ceq. So in your case, the constraint function can be defined as:
function [c, ceq] = constraint_eq(x, y, p)
% evaluate the model for required x and y
z_model = poly_model(x, y, p);
% and constrain z to be positive:
c = -z_model; % z_model >= 0, c(p) <= 0, hence c = -z_model
% no inequality constraint needed
ceq = [];
end
接下来，您需要定义一个优化函数，使实验数据与模型预测之间的误差最小化：
function err = cost_function(x, y, z, p)
z_model = poly_model(x, y, p); % determine model prediction z for x and y
ev = z_model - z; % error vector
err = norm(ev, 2)^2; % sum of squared error
end
最后，调用 fmincon 例程：
clc
clear
close all
% data
Experiment = [1.5 0.2 0.4 0.4 0.2 0.2 2.0 0.2 0.4 0.4 0.2 0.2];
Time = [0 5 10 0 5 10 0 5 10 0 5 10];
Concentration = [0 0 0 1 1 1 2 2 2 3 3 3];
% short notation for readability
x = Time;
y = Concentration;
z = Experiment;
% define XV and YV to fulfil constraint over the entire x and y domain
xv = linspace(min(x), max(x), 20);
yv = linspace(min(y), max(y), 20);
[XV, YV] = meshgrid(xv, yv);
% initial guess parameters?
p0 = rand(9, 1);
p_final = fmincon(@(p) cost_function(x, y, z, p), p0, [], [], [], [], [], [], @(p) constraint_eq(XV, YV, p));
%% check result:
ZV = poly_model(XV, YV, p_final); % evaluate points in poly23 plane
% plot result
figure(1); clf;
scatter3(x, y, z, 200, 'b.');
hold on;
surf(XV, YV, ZV)
初始参数的影响p0
正如@James Philips 在评论中指出的那样，您还可以使用无约束优化的解决方案作为约束优化的起点。对于提供的实验数据和选择的模型，你会发现实际上并没有什么区别：
% The random initial guess:
p0 = rand(9, 1);
% Optimal solution for random p0
p_rand = fmincon(@(p) cost_function(x, y, z, p), p0, [], [], [], [], [], [], @(p) constraint_eq(XV, YV, p));
% first running unconstrained optimization and use p_unc
% as start point for constrained optimization
p_unc = fmincon(@(p) cost_function(x, y, z, p), p0, [], []);
p_con= fmincon(@(p) cost_function(x, y, z, p), p_unc, [], [], [], [], [], [], @(p) constraint_eq(XV, YV, p));
% Compare errors:
SSE_unc = cost_function(x,y,z,p_unc)
SSE_con = cost_function(x,y,z,p_con)
SSE_rand = cost_function(x,y,z,p_rand)
% compare poly23 parameters
p_all = [p_unc, p_con, p_rand]
这将给出：
SSE_unc =
1.0348
SSE_con =
1.1889
SSE_rand =
1.1889
p_all =
1.3375 1.2649 1.2652
-0.3425 -0.2617 -0.2618
-1.6069 -1.0620 -1.0625
0.0258 0.0187 0.0187
0.0175 -0.0018 -0.0016
1.5708 1.0717 1.0721
-0.0042 -0.0018 -0.0018
0.0125 0.0094 0.0094
-0.3722 -0.2627 -0.2628
在这种情况下，您会发现所找到的参数存在非常小的差异，但求解器很可能需要较少的迭代次数才能获得此解。通过调整求解器设置（最优容差和约束容差），p_rand 和 p_con 的解会更接近。
一般来说，检查多个随机初始猜测是一种很好的做法，以确保您没有找到局部最小值（例如通过使用 MultiStart）。