MATLAB - NIELINIOWY ESTYMATOR LS - APROKSYMACJA

W poprzednich wpisach o aproksymacji danych pomiarowych była mowa o modelach liniowych i estymatorze LS. Często jednak jest potrzeba aproksymacji modelami nieliniowymi. W tym wpisie zajmę się tym zagadnieniem.
Model nieliniowy

We dwóch poprzednich wpisach (część 1 i część 2) na temat aproksymacji danych pomiarowych skupiłem się na modelach liniowych ze względu na współczynniki. To znaczy, że dopasowywany model do danych pomiarowych był kombinacją liniową poszukiwanych współczynników i zmiennej niezależnej x (nie występowały nieliniowości jak potęgowanie, pierwiastkowanie funkcje sinus, cosinus itd).

Było to ograniczenie wynikające z użytej postaci estymatora LS – nadaje się on wyłącznie do identyfikacji parametrów modeli liniowych.

Przykład modelu liniowego:

i nieliniowego:

Co jednak w przypadku, jeżeli dane pomiarowe muszą być opisane modelem nieliniowym ze względu na współczynniki? Na przykład funkcją wykładniczą czy sinusoidalną? W takim przypadku nie da się użyć algorytmu LS w postaci macierzowej, gdyż nie da się wprost zbudować macierzy wejść (patrz poprzedni wpis). Można jednak wykorzystać znane nam już kryterium Q (dopasowania modelu do danych pomiarowych), potraktować je jako funkcję nieznanych współczynników modelu, i numerycznie obliczyć minimum tej funkcji ze względu na te współczynniki.

Dopasowanie nieliniowego modelu, który ma aproksymować dane pomiarowe, sprowadzimy więc do problemu optymalizacji funkcji Q – znalezienia jej minimum. Wiemy przecież z poprzednich wpisów, że minimum Q, odpowiada najlepszemu dopasowaniu modelu do danych pomiarowych. Obliczone w ten sposób współczynniki modelu, które minimalizują wartość Q są poszukiwanymi przez nas współczynnikami modelu nieliniowego.

W MATLABie ten problem możemy rozwiązać na kilka sposobów, korzystając z funkcji toolboxa Optimization.

Definicja problemu

Załóżmy, że mamy dane pomiarowe o kształcie jak na rysunku. Dane można pobrać tutaj.

Od razu powinniśmy skojarzyć ten charakterystyczny kształt  - jest to krzywa Gaussa, opisana równaniem nieliniowym, o dwóch współczynnikach:

Gdzie u i s to poszukiwane przez nas współczynniki, interpretowane jako wartość oczekiwana i odchylenie standardowe. Innymi słowy u i s, to współczynniki modelu, których wartości musimy „ustalić” w taki sposób aby krzywa modelowa, jak najlepiej aproksymowała dane pomiarowe. Posłużymy się kryterium dopasowania metody LS.

Funkcja fminsearch

Pierwsze możliwe podejście to wykorzystanie funkcji fminsearch, która poszukuje minimum funkcji wielu zmiennych metodą simplexów. Nie będę wchodził w zasadę działania tej metody, zresztą nas interesuje efekt końcowy. A jest nim obliczenie takich współczynników modelu, które minimalizują Q - kryterium dopasowania modelu do danych pomiarowych.

Prześledźmy poniższy program. Dane pomiarowe można pobrać tutaj.  W pierwszej kolejności są wczytywane dane pomiarowe (wektory x i y), a następnie jest rysowany rysunek. Druga część programu to obliczenie (estymacja) poszukiwanych parametrów modelu krzywej Gaussa:

- ustalenie lub wylosowanie wartości początkowych dla u i s (jest to niezbędne w każdym zadaniu optymalizacji),

- wywołanie napisanej przez nas funkcji fitLS, w której dzieje się cała magia – o tym za chwilę,

- odebranie od funkcji fitLS wartości obliczonych współczynników modelu,

- narysowanie dopasowanej krzywej modelowej do danych pomiarowych.

Mam nadzieję, że powyższe etapy obliczeń nie wymagają wyjaśnienia. Za to takiego wyjaśnienia wymaga zawartość napisanej przez nas funkcji fitLS.

- funkcja fitLS jako argumenty przyjmuje niezbędne dane do znalezienia minimum Q: dane pomiarowe x i y, oraz wartości początkowe poszukiwanych parametrów u i s.

- funkcja fitLS zwraca nam obliczone wartości u i s zapewniające minimum kryterium Q, oraz zwraca wartość tego kryterium.

- wewnątrz fitLS w kolejnej funkcji jest zdefiniowany równaniem dopasowywany model krzywej Gaussa oraz obliczana jest wartość Q

- funkcja fminsearch tak dostraja wartości współczynników u i s aby wartość Q była minimalna.


% dane pomiarowe
load dane_fit.mat

% rysunek danych pomiarowych
figure; 
plot(x,y,'o-'); grid on;
xlabel('argument x'); ylabel('wartość y')


%% Estymacja parametrów u i s rozkładu normalnego (wg kryterium LS):

% punkt początkowy
X0 = [2 0.1];%[randn(1,1) abs(randn(1,1))];

% wywołanie funkcji fitLS, która zwraca wartości poszukiwanych parametrów
% rozkładu:
[wspl, Q] = fitLS(x, y, X0)
ue = wspl(1);
se = wspl(2);

% dopasowany model do danych pomiarowych:
yfit = (1 / (se * sqrt(2*pi))) * exp(-(x - ue).^2 ./ (2*se^2));

% rysunek
hold on; plot(x, yfit, 'LineWidth',2);
legend('dane pomiarowe', 'model')

%%

function [wspl, Q] = fitLS(x, y, X0)

    function Q = fitQ(X0)
        ue = X0(1);
        se = X0(2);

        % model 
        f = (1 / (se * sqrt(2*pi))) * exp(-(x - ue).^2 ./ (2*se^2));

%                 plot(x,y,x,f)
%                 pause(0.1)

        % kwadratowe kryterium dopasowania modelu do danych pomiarowych (LS)
        Q = sum((f-y).^2)

    end

[wspl, Q] = fminsearch(@fitQ, X0);

end


Jeżeli włączycie powyższy program, waszym oczom ukaże się rysunek jak poniżej. Poszukiwane współczynniki modelu znajdują się w zmiennej wspl.

Jeżeli jednak chcecie zobaczyć jak przebiega proces optymalizacji, i jak w kolejnych iteracjach funkcja fminsearch poszukuje najlepszego dopasowania modelu do danych, to ściągnijcie komentarz wewnątrz funkcji fitLS, tam gdzie jest plot(x,y,x,f) i pause(0.1).

Sam proces optymalizacji i działanie funkcji fminsearch może być sparametryzowane. Ale tu odsyłam Was do systemu pomocy MATLABa, gdyż dalece wykracza to poza temat tego wpisu.

Funkcja lsqcurvefit

W poprzednim punkcie sami napisaliśmy definicje kryterium Q i sami stworzyliśmy funkcje służącą do obliczania współczynników. To co ręcznie zrobiliśmy w poprzednim punkcie (dzięki temu wiemy jak to działa!) niejako automatycznie robi inna funkcja MATLABa czyli lsqcurvefit. Jej użycie jest trywialnie proste, funkcja sama oblicza wartość Q wg metody LS, a efekt jest ten sam co przy naszym poprzednim programie. Funkcja ta nie używa jednak metody simplexów, a algorytm Levenberg-Marquardt lub rust-region-reflective (do wyboru).

Użycie tej funkcji zamyka się w kilku linijkach poniżej. Myślę, że nie będziecie mieli problemu z analizą tego programu, a po tym wpisie będziecie samodzielnie dopasowywać modele nieliniowe do danych pomiarowych!


%% lsqcurvefit
clear all; close all; clc

% dane pomiarowe
load dane_fit.mat

figure; 
plot(x,y,'o-'); grid on;
xlabel('argument x'); ylabel('wartość y')

F = @(wspl,x)(1 / (wspl(2) * sqrt(2*pi))) * exp(-(x - wspl(1)).^2 ./ (2*wspl(2)^2));
X0 = [2 0.1];

[wspl_est,resnorm,~,exitflag,output] = lsqcurvefit(F,X0,x,y)
hold on
plot(x,F(wspl_est,x), 'LineWidth',2)
(Visited 482 times, 1 visits today)

Dodaj komentarz

Twój adres e-mail nie zostanie opublikowany. Wymagane pola są oznaczone *