最小二乘曲线近似

此示例演示如何使用[面向问题的优化]工作流执行非线性最小二乘曲线近似任务（https://engee.com/community/catalogs/projects/algoritm-resheniia-zadach-optimizatsii )。

我们将使用[跳转的功能。jl]图书馆(https://engee.com/helpcenter/stable/ru/julia/JuMP/index.html ）使用非线性优化库[Ipopt.jl]来制定优化问题（https://engee.com/helpcenter/stable/ru/julia/Ipopt/README.html ），随机数生成库[Random.jl]（https://engee.com/helpcenter/stable/ru/julia/stdlib/Random.html ），用于处理概率分布的库[Distributions.jl]（https://engee.com/helpcenter/stable/ru/julia/Distributions/index.html ）和库[Plots.jl]（https://engee.com/helpcenter/stable/ru/julia/juliaplots/index.html ）来可视化结果。

安装库

如果您的环境中未安装最新版本的软件包 JuMP，取消注释并运行下面的单元格:

Pkg.add(["Ipopt", "Distributions", "JuMP"])

#Pkg.add("JuMP");

要在安装完成后启动新版本的库，请单击"个人帐户"按钮:

然后点击"停止"按钮:

通过单击"Start Engee"按钮重新启动会话:

任务说明

优化任务是执行非线性最小二乘曲线近似。

此任务的模型依赖方程:

哪里 , , 和 —这些是未知参数, -响应，而t是时间。

任务的目的是找到值和其中最小化目标函数:

连接图书馆

连接库 JuMP:

using JuMP;

连接非线性求解器的库 Ipopt:

using Ipopt;

连接随机数生成库 Random:

using Random;

连接库以使用概率分布 Distributions:

using Distributions;

连接图表库 Plots:

using Plots;

数据生成

对于此任务，我们将生成人工噪声数据。

设置随机数生成器的任意粒度:

Random.seed!(0);

创建变量 A_true 和 r_true，包含模型参数的真实值。这些变量用于生成人工数据。 ydata.

对于真实值，请使用 , , 和 :

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

变量 A_true 和 r_true 它们是优化算法应该尝试重新创建的基准。将优化结果与真实值进行比较，可以评估曲线近似过程的准确性和效率。

生成200个随机值 tdata 在从0到3的范围内作为时间数据。还要创建一个变量 noisedata，包含噪声，我们将添加到最终数据 ydata.

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;

在图形上显示生成的点:

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

从图中可以看到，数据包含噪声。因此，问题的解决方案可能不会完全匹配真实参数。和 .

创建优化任务

使用函数创建优化任务 Model() 并在括号中指定求解器的名称。:

model = Model(Ipopt.Optimizer)

A JuMP Model
├ solver: Ipopt
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
└ Names registered in the model: none

创建变量和每个包含两个值, – , , 和 :

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

创建一个最小化平方和的非线性目标函数:

@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)));

设置适当的初始值是优化问题高效准确解决的重要步骤，显着影响优化过程的速度、质量和可靠性。

指定初始值 , , 和 .

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

打印初始值和 :

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

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

解决问题

解决优化问题:

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_sol = value.(A);
r_sol = value.(r);

输出优化结果:

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

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

结果的可视化和分析

计算曲线近似的结果:

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

在图形上显示结果:

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

您可以将获得的结果与最初设置的值进行比较。 , , 和 . 并且round函数将允许您四舍五入过度准确的值。

比较结果与原始值:

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

比较结果与原始值:

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

所以结果是平均而言，它们与初始值的不同之处在于，而结果上 .

结论

在这个例子中，我们使用面向问题的方法使用最小二乘法（OLS）解决了曲线近似的非线性问题。我们还将结果与模型的真实参数进行了可视化和比较。