Нелинейный метод наименьших квадратов, основанный на проблеме

В этом примере показано, как выполнить аппроксимирование кривыми нелинейного метода наименьших квадратов с помощью Основанного на проблеме Рабочего процесса Оптимизации.

Модель

Уравнение модели для этой проблемы

$y (t) = A_{1} \exp (r_{1} t) + A_{2} \exp (r_{2} t),$

где $A_{1}$ , $A_{2}$ , $r_{1}$ , и $r_{2}$ неизвестные параметры, $y$ ответ, и $t$ время. Проблема требует данных в течение многих времен tdata и (шумные) измерения ответа ydata. Цель состоит в том, чтобы найти лучшее $A$ и $r$ , значение тех значений, которые минимизируют

$\sum_{t \in t d a t a} {(y (t) - y d a t a)}^{2} .$

Выборочные данные

Как правило, у вас есть данные для проблемы. В этом случае сгенерируйте искусственные зашумленные данные для проблемы. Используйте 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);

Создайте выражение для целевой функции, которая является суммой квадратов, чтобы минимизировать. Поскольку ответ имеет экспоненциальный термин, запишите анонимную функцию для ответа и преобразуйте ответ на выражение оптимизации с помощью 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);

Создайте задачу оптимизации с целевой функцией 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((anonymousFunction1(A, r) - extraParams{1}).^2)

       where:

         anonymousFunction1 = @(A,r)A(1)*exp(r(1)*tdata)+A(2)*exp(r(2)*tdata);

         extraParams{1}:

         2.9278
         2.7513
         2.7272
         2.4199
         2.3172
         2.3961
         2.2522
         2.1974
         2.1666
         2.0944
         1.9566
         1.7989
         1.7984
         1.7540
         1.8318
         1.6745
         1.6874
         1.5526
         1.5229
         1.5680
         1.4784
         1.5280
         1.3727
         1.2968
         1.4012
         1.3602
         1.2714
         1.1773
         1.2119
         1.2033
         1.2037
         1.1729
         1.1829
         1.1602
         1.0448
         1.0320
         1.0397
         1.0334
         1.0233
         1.0275
         0.8173
         0.9373
         1.0202
         0.8896
         0.9791
         0.9128
         0.7763
         0.7669
         0.6579
         0.7135
         0.7978
         0.7164
         0.7071
         0.6429
         0.6676
         0.6782
         0.6802
         0.6328
         0.6301
         0.7406
         0.4908
         0.7126
         0.6136
         0.6269
         0.4668
         0.4963
         0.5007
         0.5226
         0.3764
         0.4824
         0.3930
         0.4390
         0.4665
         0.4490
         0.4841
         0.4539
         0.3698
         0.3974
         0.3356
         0.3045
         0.4131
         0.3561
         0.3506
         0.3960
         0.3625
         0.3446
         0.3778
         0.3565
         0.3187
         0.2677
         0.2664
         0.3572
         0.2129
         0.2919
         0.1764
         0.3210
         0.3016
         0.2572
         0.2514
         0.1301
         0.2825
         0.1372
         0.1243
         0.2421
         0.1888
         0.2547
         0.2559
         0.2632
         0.1801
         0.2309
         0.2134
         0.2495
         0.2332
         0.2512
         0.1875
         0.1861
         0.2397
         0.0803
         0.1579
         0.1196
         0.1541
         0.1978
         0.2034
         0.1095
         0.1332
         0.1567
         0.1345
         0.1635
         0.1661
         0.0991
         0.1366
         0.0387
         0.1922
         0.1031
         0.0714
         0.1178
         0.0568
         0.1255
         0.0957
         0.2313
         0.1710
        -0.0148
         0.1316
         0.0385
         0.0946
         0.1147
         0.1436
         0.0917
         0.1840
         0.0786
         0.1161
         0.1327
         0.1026
         0.1421
         0.1142
         0.0553
         0.0036
         0.1866
         0.0634
         0.0974
         0.1203
         0.0939
         0.0429
         0.0640
         0.0811
         0.1603
         0.0427
         0.1244
         0.0993
         0.0696
         0.0264
         0.0641
         0.0703
         0.0010
         0.0793
         0.0267
         0.0625
         0.0834
         0.0204
         0.0507
         0.0826
        -0.0272
         0.1161
         0.1832
         0.1100
         0.0453
         0.0826
         0.0079
         0.1531
         0.1052
         0.0965
         0.0132
         0.0729
         0.0287
         0.0410
         0.0280
         0.0049
         0.0102
         0.0442
        -0.0343

Основанное на проблеме решение

Решите задачу.

[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.

sol = struct with fields:
    A: [2x1 double]
    r: [2x1 double]

fval = 0.4724

Постройте получившееся решение и исходные данные.

figure
responsedata = evaluate(response,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.

Смотрите также

solve

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