Нелинейный метод наименьших квадратов, основанный на проблеме

В этом примере показано, как выполнить аппроксимирование кривыми нелинейного метода наименьших квадратов с помощью Основанного на проблеме Рабочего процесса Оптимизации.

Модель

Уравнение модели для этой проблемы

y(t)=A1exp(r1t)+A2exp(r2t),

где A1, A2, r1, и r2 неизвестные параметры, y ответ, и t время. Проблема требует данных в течение многих времен tdata и (шумные) измерения ответа ydata. Цель состоит в том, чтобы найти лучшее A и r, значение тех значений, которые минимизируют

ttdata(y(t)-ydata)2.

Выборочные данные

Как правило, у вас есть данные для проблемы. В этом случае сгенерируйте искусственные зашумленные данные для проблемы. Используйте A = [1,2] и r = [-1,-3] как базовые значения и использование 200 случайных значений от 0 к 3 как данные времени. Постройте получившиеся точки данных.

rng default % For reproducibility
A = [1,2];
r = [-1,-3];
tdata = 3*rand(200,1);
tdata = sort(tdata); % Increasing times for easier plotting
noisedata = 0.05*randn(size(tdata)); % Artificial noise
ydata = A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata) + noisedata;
plot(tdata,ydata,'r*')
xlabel 't'
ylabel 'Response'

Данные являются шумными. Поэтому решение, вероятно, не будет совпадать с исходными параметрами A и r очень хорошо.

Подход, основанный на проблеме

Найти параметры оптимальной подгонки A и r, сначала задайте переменные оптимизации с теми именами.

A = optimvar('A',2);
r = optimvar('r',2);

Создайте выражение для целевой функции, которая является суммой квадратов, чтобы минимизировать.

fun = A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata);
obj = sum((fun - ydata).^2);

Создайте задачу оптимизации с целевой функцией obj.

lsqproblem = optimproblem("Objective",obj);

Для подхода, основанного на проблеме задайте начальную точку как структуру с именами переменных как поля структуры. Задайте начальный A = [1/2,3/2] и начальный r = [-1/2,-3/2].

x0.A = [1/2,3/2];
x0.r = [-1/2,-3/2];

Рассмотрите формулировку задачи.

show(lsqproblem)
  OptimizationProblem : 

	Solve for:
       A, r

	minimize :
       sum(arg6)

       where:

         arg5 = extraParams{3};
         arg6 = (((A(1) .* exp((r(1) .* extraParams{1}))) + (A(2) .* exp((r(2) .* extraParams{2})))) - arg5).^2;

       extraParams

Основанное на проблеме решение

Решите задачу.

[sol,fval] = solve(lsqproblem,x0)
Solving problem using lsqnonlin.

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

<stopping criteria details>
sol = struct with fields:
    A: [2×1 double]
    r: [2×1 double]

fval = 0.4724

Постройте получившееся решение и исходные данные.

figure
responsedata = evaluate(fun,sol);
plot(tdata,ydata,'r*',tdata,responsedata,'b-')
legend('Original Data','Fitted Curve')
xlabel 't'
ylabel 'Response'
title("Fitted Response")

График показывает, что подходящие данные совпадают с исходными зашумленными данными довольно хорошо.

Смотрите, как тесно подходящие параметры совпадают с исходными параметрами A = [1,2] и r = [-1,-3].

disp(sol.A)
    1.1615
    1.8629
disp(sol.r)
   -1.0882
   -3.2256

Подходящие параметры выключены приблизительно на 15% в A и 8% в r.

Неподдерживаемые функции требуют fcn2optimexpr

Если ваша целевая функция не состоит из элементарных функций, необходимо преобразовать функцию в выражение оптимизации с помощью fcn2optimexpr. Смотрите Преобразуют Нелинейную Функцию в Выражение Оптимизации. Для существующего примера:

fun = @(A,r) A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata);
response = fcn2optimexpr(fun,A,r);
obj = sum((response - ydata).^2);

Остаток от шагов в решении задачи является тем же самым. Единственное другое различие находится в стандартной программе графического вывода, где вы вызываете response вместо fun:

responsedata = evaluate(response,sol);

Для списка поддерживаемых функций смотрите Поддерживаемые Операции на Переменных и выражениях Оптимизации.

Смотрите также

Похожие темы