Engee 文档
Notebook

使用最小二乘法进行曲线逼近

本例演示了如何使用面向问题的优化 工作流程执行非线性曲线逼近（近似）问题。

我们将使用JuMP.jl库的函数来提出优化问题，非线性优化库Ipopt.jl、随机数生成库Random.jl、用于处理概率分布的库Distributions.jl和库Plots.jl来可视化结果。

安装库

如果您的环境中没有安装最新版本的JuMP 软件包，请取消注释并运行下面的方框：

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

要在安装完成后启动新版本的程序库，请单击 "我的账户 "按钮：

screenshot_20240710_134852.png 然后点击 "停止 "按钮：

screenshot_20240710_2.png

按下 "Start 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 个从 0 到 3 的tdata 随机值作为时间数据。同时创建一个包含噪声的变量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$ 。四舍五入功能可以对过于精确的数值进行四舍五入。

将$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\%$ 。

结论

在本例中，我们采用以问题为导向的方法，使用最小二乘法 (LSM) 解决了一个非线性曲线逼近问题。我们还将结果可视化，并与真实的模型参数进行了比较。