MATLAB Curve Fitting в 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 давать мне приближение, которое всегда дает Z> 0, если времяот 0 до 10, а концентрация от 0 до 3?

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) - массив нелинейных ограничений неравенства приИкс. fmincon пытается удовлетворить

  c (x) <= 0 для всех элементов c. </p>

• ceq (x) - это массив нелинейных ограничений равенства в точке x. fmincon пытается удовлетворить

  ceq (x) = 0 для всех записей ceq. Таким образом, в вашем случае функцию ограничения можно определить следующим образом:

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).