MATLAB 3D曲线拟合，带有附加边界

Question

简介

假设我有一组实验数据，需要找到一个描述所选序列的多项式逼近。实验结果取决于两个变量-时间和浓度。让示例数据如下所示：

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时，多项式结果为负。我正在研究的系统是生物系统。因此，从物理上讲不可能有任何负值。所以我的问题是：如何设置边界/约束以防止所得多项式在实验范围内为负？如何强制MATLAB给我近似值，如果时间在0与0之间，则Z总是为0。和10，浓度在0和3之间？

Answer 1

对于使用fmincon的非线性约束优化，首先需要定义一个确定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的矩阵输入函数，这是计算要在实验数据范围内施加的约束所必需的。

约束已在单独的函数中定义。从文档：

非线性约束，指定为函数句柄或函数名称。 nonlcon是一个接受向量x或数组x并返回两个数组c（x）和ceq（x）的函数。

• c（x）是x处的非线性不等式约束的数组。 fmincon试图满足

  c（x）<= 0对于c的所有条目。

• ceq（x）是x处的非线性等式约束的数组。 fmincon试图满足

  ceq（x）= 0对于ceq的所有条目。因此，在您的情况下，约束函数可以定义为：

y

接下来，您需要定义一个优化函数，以最小化实验数据与模型预测之间的误差：

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)

初始参数的影响

正如@James Philips在评论中指出的那样，您还可以将无约束优化的解决方案用作约束优化的起点。对于提供的实验数据和选择的模型，您将看到并没有真正的区别：

p0

哪个会给：

% 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）。