Аппроксимация нелинейной функции
artpgchНелинейная аппроксимация данных
Введение
В данном примере мы рассмотрим, как аппроксимировать нелинейную функцию по данным с использованием различных алгоритмов Julia.
Присоединим необходимые библиотеки:
using Plots, LsqFit, Optim
Постановка задачи
Задача состоит в аппроксимации функции по следующим данным:
xdata = [
0.0000 5.8955
0.1000 3.5639
0.2000 2.5173
0.3000 1.9790
0.4000 1.8990
0.5000 1.3938
0.6000 1.1359
0.7000 1.0096
0.8000 1.0343
0.9000 0.8435
1.0000 0.6856
1.1000 0.6100
1.2000 0.5392
1.3000 0.3946
1.4000 0.3903
1.5000 0.5474
1.6000 0.3459
1.7000 0.1370
1.8000 0.2211
1.9000 0.1704
2.0000 0.2636
];
Построим график точек данных:
t = xdata[:, 1];
y = xdata[:, 2];
график1 = plot();
scatter!(t, y, title = "Точки данных", markerstrokecolor = :blue, markercolor = :white, markerstrokewidth = 1, markersize = 4, legend = false);
display(график1)
Попробуем аппроксимировать функцию к данным на графике.
Решение с помощью функции curve_fit
Чтобы аппроксимировать функцию к данным с помощью curve_fit, сопоставим параметры в функции с переменной x следующим образом:
x(1) =
x(2) =
x(3) =
x(4) =
Определим функцию для аппроксимации:
F(xdata, x) = x[1] * exp.(-x[2] .* xdata) + x[3] * exp.(-x[4] .* xdata);
Произвольно зададим начальную точку x0 следующим образом:
x0 = [1.0, 1.0, 1.0, 0.0];
Запустим функцию curve_fit и построим график полученной аппроксимации:
start_time = time()
fit = curve_fit(F, t, y, x0);
x = fit.param;
resnorm = sum(fit.resid.^2);
elapsed_time = time() - start_time;
println(x)
println("Время вычисления: $elapsed_time секунд")
график2 = scatter(t, y, markerstrokecolor=:blue, markercolor=:white, markerstrokewidth=1, markersize=4, label="Точки данных")
plot!(t, F(t, x), linewidth=2, color=:blue, label="Аппроксимация")
title!("Аппроксимация")
display(график2)
Решение с помощью функции optimize
Чтобы решить задачу с помощью optimize, определим целевую функцию как сумму квадратов остатков:
Fsumsquares(x) = sum((F(t, x) - y).^2);
start_time = time()
result = optimize(Fsumsquares, x0, LBFGS());
xunc = Optim.minimizer(result) ;
ressquared = Optim.minimum(result) ;
eflag = Optim.converged(result) ? 1 : 0 ;
elapsed_time = time() - start_time;
println(xunc)
println("Время вычисления: $elapsed_time секунд")
Заметим, что optimize находит то же решение, что и curve_fit, но требует значительно больше времени для вычислений. Так как порядок переменных произвольный, результаты вычислений расположены в ином порядке.
Разделение линейной и нелинейной задач
Заметим, что задача аппроксимации линейна по параметрам и . Это означает, что для любых значений и можно использовать оператор обратного деления для нахождения значений c₁ и c₂, решающих задачу наименьших квадратов.
Переформулируем задачу как двумерную, выполняя поиск оптимальных значений и . Значения и вычисляются на каждом шаге с использованием оператора обратного деления.
function fitvector(lam, xdata, ydata)
A = [exp(-λ * x) for x in xdata, λ in lam]
c = A \ ydata
yEst = A * c
return yEst
end
Решим задачу с использованием curve_fit, начиная с двумерной начальной точки [1, 0]:
x02 = [1.0, 0.0];
F2(t_data, x) = fitvector(x, t_data, y)
start_time = time()
fit2 = curve_fit(F2, t, y, x02);
x2 = fit2.param;
resnorm2 = sum(fit2.resid.^2);
exitflag2 = fit2.converged ? 1 : 0;
elapsed_time = time() - start_time;
println(x2)
#println("Время вычисления: $elapsed_time секунд")
Разделённая задача наиболее устойчива к начальному приближению
Выбор неудачной начальной точки для исходной задачи с четырьмя параметрами приводит к локальному решению, которое не является глобальным. Однако выбор начальной точки с такими же неудачными значениями λ₁ и λ₂ для разделённой задачи с двумя параметрами, приводит к глобальному решению. Чтобы продемонстрировать этот результат, перезапустим исходную задачу с начальной точкой, которая приводит к относительно неудачному локальному решению, и сравним полученную аппроксимацию с глобальным решением.
x0bad = [10.0, 1.0, 1.0, 0.0]
start_time = time()
fitbad = curve_fit(F, t, y, x0bad);
xbad = fitbad.param;
resnormbad = sum(fitbad.resid.^2);
exitflagbad = fitbad.converged ? 1 : 0;
elapsed_time = time() - start_time;
println(xbad)
println("Время вычисления: $elapsed_time секунд")
график3 = scatter(t, y, markerstrokecolor=:blue, markercolor=:white, markerstrokewidth=1, markersize=4, label="Данные")
plot!(t, F(t, x), linewidth=2, color=:blue, label="Удачные параметры")
plot!(t, F(t, xbad), linewidth=2, color=:red, label="Неудачные параметры")
title!("Аппроксимация")
display(график3)
println("Остаточная норма в удачной конечной точке: $(resnorm);")
println("Остаточная норма в неудачной конечной точке $(resnormbad).")
Вывод
Данный пример демонстрирует ключевые аспекты нелинейной аппроксимации данных с использованием методов оптимизации:
-
Эффективность специализированных алгоритмов - функция curve_fit, предназначенная специально для задач аппроксимации кривых, показала значительно более высокую эффективность по сравнению с функцией optimize, требуя меньшего количества вычислений целевой функции для достижения того же результата.
-
Важность разделения параметров - выделение линейных параметров задачи позволяет существенно упростить исходную четырехпараметрическую задачу до двумерной, что не только сохраняет эффективность вычислений, но и повышает устойчивость метода.
-
Устойчивость к начальным приближениям - разделённая формулировка задачи продемонстрировала повышенную надёжность, позволяя избегать локальных минимумов и находить глобальное решение даже при неудачном выборе начальных значений параметров.
-
Практическая значимость - представленные подходы имеют широкую применимость в различных областях, где требуется точная аппроксимация экспериментальных данных сложными нелинейными моделями, включая химическую кинетику, биологические процессы и физические эксперименты.
Метод разделения линейных и нелинейных параметров представляет особую ценность, поскольку сочетает вычислительную эффективность с повышенной надёжностью сходимости к глобальному решению, что делает его предпочтительным выбором для решения практических задач нелинейной аппроксимации.