Нелинейный метод наименьших квадратов, основанный на проблеме
В этом примере показано, как выполнить аппроксимирование кривыми нелинейного метода наименьших квадратов с помощью Основанного на проблеме Рабочего процесса Оптимизации.
Модель
Уравнение модели для этой проблемы
где , , , и неизвестные параметры, ответ, и время. Проблема требует данных в течение многих времен tdata и (шумные) измерения ответа ydata. Цель состоит в том, чтобы найти лучшее и , значение тех значений, которые минимизируют
Выборочные данные
Как правило, у вас есть данные для проблемы. В этом случае сгенерируйте искусственные зашумленные данные для проблемы. Используйте 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);
Создайте выражение для целевой функции, которая является суммой квадратов, чтобы минимизировать.
fun = A(1)*exp(r(1)*tdata) + A(2)*exp(r(2)*tdata); obj = sum((fun - 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(arg6)
where:
arg5 = extraParams{3};
arg6 = (((A(1) .* exp((r(1) .* extraParams{1})))
+ (A(2) .* exp((r(2) .* extraParams{2})))) - arg5).^2;
extraParams{1}:
0.0139
0.0357
0.0462
0.0955
0.1033
0.1071
0.1291
0.1385
0.1490
0.1619
0.1793
0.2276
0.2279
0.2345
0.2434
0.2515
0.2533
0.2894
0.2914
0.2926
0.3200
0.3336
0.3570
0.3700
0.3810
0.3897
0.3959
0.4082
0.4159
0.4257
0.4349
0.4366
0.4479
0.4571
0.4728
0.4865
0.4878
0.4969
0.5070
0.5136
0.5455
0.5505
0.5517
0.5606
0.5669
0.5898
0.6714
0.6869
0.7043
0.7197
0.7199
0.7251
0.7306
0.7533
0.7628
0.7653
0.7725
0.7796
0.7889
0.7914
0.8281
0.8308
0.8355
0.8575
0.8890
0.9190
0.9336
0.9513
1.0114
1.0132
1.0212
1.0500
1.0529
1.0550
1.0595
1.1055
1.1077
1.1413
1.1447
1.1692
1.1767
1.1993
1.2054
1.2117
1.2518
1.2653
1.2942
1.3076
1.3162
1.3280
1.3368
1.3404
1.3516
1.3528
1.4082
1.4199
1.4561
1.4604
1.4678
1.4693
1.4726
1.4951
1.5179
1.5255
1.5323
1.5397
1.5856
1.5924
1.6150
1.6406
1.6410
1.6416
1.6492
1.6496
1.7035
1.7065
1.7256
1.7391
1.7558
1.7558
1.8059
1.8481
1.8662
1.8769
1.8971
1.9389
1.9432
1.9473
1.9622
1.9653
1.9664
1.9672
2.0362
2.0391
2.0603
2.0676
2.0845
2.0972
2.1181
2.1281
2.1952
2.2294
2.2341
2.2445
2.2538
2.2612
2.2641
2.2716
2.2732
2.2966
2.3247
2.3271
2.3375
2.3407
2.3408
2.3766
2.3829
2.3845
2.3856
2.4002
2.4008
2.4429
2.4442
2.4519
2.4529
2.4636
2.4704
2.4775
2.4925
2.5222
2.5474
2.5591
2.6061
2.6079
2.6727
2.7002
2.7081
2.7174
2.7319
2.7400
2.7401
2.7472
2.7516
2.7878
2.7882
2.8020
2.8020
2.8262
2.8344
2.8507
2.8684
2.8715
2.8725
2.8779
2.8785
2.8792
2.8857
2.8947
2.9118
2.9884
extraParams{2}:
0.0139
0.0357
0.0462
0.0955
0.1033
0.1071
0.1291
0.1385
0.1490
0.1619
0.1793
0.2276
0.2279
0.2345
0.2434
0.2515
0.2533
0.2894
0.2914
0.2926
0.3200
0.3336
0.3570
0.3700
0.3810
0.3897
0.3959
0.4082
0.4159
0.4257
0.4349
0.4366
0.4479
0.4571
0.4728
0.4865
0.4878
0.4969
0.5070
0.5136
0.5455
0.5505
0.5517
0.5606
0.5669
0.5898
0.6714
0.6869
0.7043
0.7197
0.7199
0.7251
0.7306
0.7533
0.7628
0.7653
0.7725
0.7796
0.7889
0.7914
0.8281
0.8308
0.8355
0.8575
0.8890
0.9190
0.9336
0.9513
1.0114
1.0132
1.0212
1.0500
1.0529
1.0550
1.0595
1.1055
1.1077
1.1413
1.1447
1.1692
1.1767
1.1993
1.2054
1.2117
1.2518
1.2653
1.2942
1.3076
1.3162
1.3280
1.3368
1.3404
1.3516
1.3528
1.4082
1.4199
1.4561
1.4604
1.4678
1.4693
1.4726
1.4951
1.5179
1.5255
1.5323
1.5397
1.5856
1.5924
1.6150
1.6406
1.6410
1.6416
1.6492
1.6496
1.7035
1.7065
1.7256
1.7391
1.7558
1.7558
1.8059
1.8481
1.8662
1.8769
1.8971
1.9389
1.9432
1.9473
1.9622
1.9653
1.9664
1.9672
2.0362
2.0391
2.0603
2.0676
2.0845
2.0972
2.1181
2.1281
2.1952
2.2294
2.2341
2.2445
2.2538
2.2612
2.2641
2.2716
2.2732
2.2966
2.3247
2.3271
2.3375
2.3407
2.3408
2.3766
2.3829
2.3845
2.3856
2.4002
2.4008
2.4429
2.4442
2.4519
2.4529
2.4636
2.4704
2.4775
2.4925
2.5222
2.5474
2.5591
2.6061
2.6079
2.6727
2.7002
2.7081
2.7174
2.7319
2.7400
2.7401
2.7472
2.7516
2.7878
2.7882
2.8020
2.8020
2.8262
2.8344
2.8507
2.8684
2.8715
2.8725
2.8779
2.8785
2.8792
2.8857
2.8947
2.9118
2.9884
extraParams{3}:
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(fun,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.
Неподдерживаемые функции требуют fcn2optimexpr
Если ваша целевая функция не состоит из элементарных функций, необходимо преобразовать функцию в выражение оптимизации с помощью 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);
Остаток от шагов в решении задачи является тем же самым. Единственное другое различие находится в стандартной программе графического вывода, где вы вызываете response вместо fun:
responsedata = evaluate(response,sol);
Для списка поддерживаемых функций смотрите Поддерживаемые Операции для Переменных и выражений Оптимизации.