В этом примере показов, как выполнить нелинейный метод наименьших квадратов аппроксимирования кривыми с помощью рабочего процесса оптимизации на основе проблем.
Модельное уравнение для этой задачи
где , , , и являются неизвестными параметрами, является ответом, и время. Задача требует данных для времени tdata
и (шумные) измерения отклика ydata
. Цель - найти лучшее и , что означает те значения, которые минимизируют
Как правило, у вас есть данные для задачи. В этом случае сгенерируйте искусственные зашумленные данные для задачи. Использование 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);
Список поддерживаемых функций см. в Поддерживаемые операции с переменными оптимизации и выражениями.