带有附加边界的三维MATLAB曲线拟合

Question

Introduction

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

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到10之间，浓度在0到3之间时，它总是给出Z>0？

Answer 1

对于使用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的矩阵输入函数，这是计算要在实验数据范围内施加的约束所必需的。

约束已在单独的函数中定义。从医生那里：

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

• c(x)是在x.fmincon试图满足的非线性不等式约束的数组。

c(x) <= 0对于c.

• ceq(x)的所有条目都是在x.fmincon试图满足的非线性等式约束的数组。

对于ceq的所有条目，ceq(x) =0。因此，在您的示例中，约束函数可以定义为：

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在评论中指出的那样，您也可以使用无约束优化的解决方案作为约束优化的起点。对于所提供的实验数据和选择的模型，您将看到并没有真正的区别：

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