Документация Engee
Notebook

Приближение кривой методом наименьших квадратов

Открыть пример в Engee

В данном примере продемонстрировано, как выполнить нелинейную задачу по приближению (аппроксимации) кривой методом наименьших квадратов (МНК) с помощью рабочего процесса проблемно-ориентированной оптимизации.

Мы воспользуемся функциями библиотеки JuMP.jl для формулировки оптимизационной задачи, библиотекой нелинейной оптимизации Ipopt.jl, библиотекой генерации случайных чисел Random.jl, библиотекой для работы с распределениями вероятностей Distributions.jl и библиотекой Plots.jl для визуализации результатов.

Установка библиотек

Если в вашем окружении не установлена последняя версия пакета JuMP, раскомментируйте и запустите ячейку ниже:

In [ ]: 
#Pkg.add("JuMP");

Для запуска новой версии библиотеки после завершения установки нажмите на кнопку «Личный кабинет»:

screenshot_20240710_134852.png Затем нажмите на кнопку «Стоп»:

screenshot_20240710_2.png

Перезапустите сессию нажатием кнопки «Старт Engee»:

screenshot_20240710_135431.png

Описание задачи

Задача оптимизации состоит в том, чтобы выполнить нелинейную аппроксимацию кривой методом наименьших квадратов.

Модельное уравнение зависимости для этой задачи:

$$ y(t)=A_1exp(r_1t)+A_2exp(r_2t) $$

где $A_1$, $A_2$, $r_1$ и $r_2$ — это неизвестные параметры,$y$ — отклик, а t — время.

Цель задачи состоит в том, чтобы найти значения $A$ и $r$, которые минимизируют целевую функцию:

$$ \sum \limits _{t \in \mathbb{tdata} } (y(t)-ydata)^2 $$

Подключение библиотек

Подключите библиотеку JuMP:

In [ ]: 
using JuMP;

Подключите библиотеку нелинейного решателя Ipopt:

In [ ]: 
using Ipopt;

Подключите библиотеку генерации случайных чисел Random:

In [ ]: 
using Random;

Подключите библиотеку для работы с распределениями вероятностей Distributions:

In [ ]: 
using Distributions;

Подключите библиотеку построения графиков Plots:

In [ ]: 
using Plots;

Генерация данных

Для данной задачи мы сгенерируем искусственные зашумленные данные.

Задайте произвольное зерно генератора случайных чисел:

In [ ]: 
Random.seed!(0);

Создайте переменные A_true и r_true, содержащие истинные значения параметров модели. Эти переменные используются для генерации искусственных данных ydata.

В качестве истинных значений используйте $A_1=1$, $A_2=2$, $r_1=-1$ и $r_2=-3$:

In [ ]: 
A_true = [1, 2];
r_true = [-1, -3];

Переменные A_true и r_true являются эталоном, который алгоритм оптимизации должен попытаться воссоздать. Сравнение результатов оптимизации с истинными значениями позволяет нам оценить точность и эффективность процесса аппроксимации кривой.

Сгенерируйте 200 случайных значений tdata в диапазоне от 0 до 3 в качестве данных времени. Также создайте переменную noisedata, содержащую шум, которую мы добавим к итоговым данным ydata.

In [ ]: 
tdata = sort(3 * rand(200));
noisedata = 0.05 * randn(length(tdata)); 
ydata = A_true[1] * exp.(r_true[1] * tdata) + A_true[2] * exp.(r_true[2] * tdata) + noisedata;

Выведите полученные точки на график:

In [ ]: 
plot(tdata, ydata, seriestype=:scatter, color=:red, legend=:right, label="Исходные данные")
xlabel!("t")
ylabel!("Отклик")
Out[0]:

Как видно из графика, данные содержат шум. Поэтому решение задачи, вероятно, не будет точно совпадать с истинными параметрами $A$ и $r$.

Создание задачи оптимизации

Создайте оптимизационную задачу с помощью функции Model() и укажите название решателя в скобках:

In [ ]: 
model = Model(Ipopt.Optimizer)
Out[0]:
A JuMP Model
├ solver: Ipopt
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

Создайте переменные $A$ и $r$, содержащие по два значения, – $A_1$, $A_2$, $r_1$ и $r_2$:

In [ ]: 
@variable(model, A[1:2]);
@variable(model, r[1:2]);

Создайте нелинейную целевую функцию, минимизирующую сумму квадратов:

In [ ]: 
@NLobjective(model, Min, sum((A[1]*exp(r[1]*t) + A[2]*exp(r[2]*t) - y)^2 for (t,y) in zip(tdata, ydata)));

Установка соответствующих начальных значений – важный шаг для эффективного и точного решения задач оптимизации, существенно влияющий на скорость, качество и надежность процесса оптимизации.

Укажите начальные значения $A_1=\frac{1}{2}$, $A_2=\frac{3}{2}$, $r_1=-\frac{1}{2}$ и $r_2=-\frac{3}{2}$.

In [ ]: 
set_start_value.(A, [1/2, 3/2]);
set_start_value.(r, [-1/2, -3/2]);

Выведите начальные значения $A$ и $r$:

In [ ]: 
println("Начальные значения A: ", start_value.(A))
println("Начальные значения r: ", start_value.(r))

Начальные значения A: [0.5, 1.5]
Начальные значения r: [-0.5, -1.5]

Решение задачи

Решите оптимизационную задачу:

In [ ]: 
optimize!(model)

This is Ipopt version 3.14.13, running with linear solver MUMPS 5.6.1.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       10

Total number of variables............................:        4
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  7.7600291e+00 0.00e+00 1.27e+01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  5.9773864e+00 0.00e+00 8.83e+00  -1.0 8.14e-02   2.0 1.00e+00 1.00e+00f  1
   2  3.5889330e+00 0.00e+00 5.31e+00  -1.0 1.50e-01   1.5 1.00e+00 1.00e+00f  1
   3  1.6781193e+00 0.00e+00 2.57e+00  -1.0 1.98e-01   1.0 1.00e+00 1.00e+00f  1
   4  8.3132157e-01 0.00e+00 1.17e+00  -1.0 2.34e-01   0.6 1.00e+00 1.00e+00f  1
   5  5.4494882e-01 0.00e+00 4.52e-01  -1.0 2.47e-01   0.1 1.00e+00 1.00e+00f  1
   6  4.9173214e-01 0.00e+00 5.15e-01  -1.7 2.48e-01    -  1.00e+00 1.00e+00f  1
   7  4.8601288e-01 0.00e+00 3.37e-02  -1.7 5.88e-02  -0.4 1.00e+00 1.00e+00f  1
   8  4.8080205e-01 0.00e+00 3.94e-01  -2.5 2.66e-01    -  1.00e+00 1.00e+00f  1
   9  4.7632633e-01 0.00e+00 8.84e-02  -2.5 1.17e-01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  4.7615210e-01 0.00e+00 3.43e-03  -2.5 3.40e-02    -  1.00e+00 1.00e+00f  1
  11  4.7615108e-01 0.00e+00 1.79e-05  -3.8 2.09e-03    -  1.00e+00 1.00e+00f  1
  12  4.7615108e-01 0.00e+00 2.24e-09  -8.6 3.05e-05    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 12

                                   (scaled)                 (unscaled)
Objective...............:   4.7615108059062977e-01    4.7615108059062977e-01
Dual infeasibility......:   2.2379096544650201e-09    2.2379096544650201e-09
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.2379096544650201e-09    2.2379096544650201e-09


Number of objective function evaluations             = 13
Number of objective gradient evaluations             = 13
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 12
Total seconds in IPOPT                               = 0.010

EXIT: Optimal Solution Found.

Сохраните значения $A$ и $r$ в переменных:

In [ ]: 
A_sol = value.(A);
r_sol = value.(r);

Выведите результаты оптимизации:

In [ ]: 
println("Оптимизированные значения A: ", A_sol)
println("Оптимизированные значения r: ", r_sol)

Оптимизированные значения A: [1.0151611001028849, 2.0255646656856654]
Оптимизированные значения r: [-1.0125976166336879, -3.0589072699017805]

Визуализация и анализ результатов

Рассчитайте результат аппроксимации кривой:

In [ ]: 
fitted_response = A_sol[1] * exp.(r_sol[1] * tdata) + A_sol[2] * exp.(r_sol[2] * tdata);

Выведите результаты на график:

In [ ]: 
plot!(tdata, fitted_response, color=:blue, label="Приближённая кривая")
title!("Приближённый отклик")
Out[0]:

Вы можете сравнить полученные результаты с изначально заданными значения $A_1=1$, $A_2=2$, $r_1=-1$ и $r_2=-3$. А функция round позволит округлить избыточно точные значения.

Сравните результаты $A$ с изначальными значениями:

In [ ]: 
percent_diffs = ((A_sol .- A_true) ./ A_true) .* 100
avg_percent_diff = mean(percent_diffs)
println("Процентная разница: ",  round.(percent_diffs, digits=2))
println("Средняя процентная разница: ", round.(avg_percent_diff, digits=2))

Процентная разница: [1.52, 1.28]
Средняя процентная разница: 1.4

Сравните результаты $r$ с изначальными значениями:

In [ ]: 
percent_diffs = ((r_sol .- r_true) ./ r_true) .* 100
avg_percent_diff = mean(percent_diffs)
println("Процентная разница: ", round.(percent_diffs, digits=2))
println("Средняя процентная разница: ", round.(avg_percent_diff, digits=2))

Процентная разница: [1.26, 1.96]
Средняя процентная разница: 1.61

Таким образом, результаты $A$ в среднем отличаются от изначальных значений на $1.4\%$, а результаты $r$ на $1.61\%$.

Заключение

В данном примере мы решили нелинейную задачу по приближению (аппроксимации) кривой с помощью метода наименьших квадратов (МНК), используя проблемно-ориентированный подход. Мы также визуализировали и произвели сравнение результатов с истинными параметрами модели.