Moduł 12 · identyfikacja

Identyfikacja parametryczna

Fitowanie modelu do pomiarów, ARX/ARMAX, RLS online, persistent excitation, closed-loop ID, praktyczny workflow.

1. Po co identyfikacja — most teoria↔praktyka

Analogia. Dostajesz na biurko nowy instrument — mówią Ci tylko, że to gitara. Nie znasz napięcia strun, długości pudła, materiału. Nie zagrasz gitarą, nie wiedząc tego. Identyfikacja parametryczna to strojenie instrumentu przed użyciem: wzbudzasz znanym sygnałem (palcami trącasz strunę), mierzysz odpowiedź (słuchasz tonu), wyciągasz parametry (długość, napięcie).

Wszystkie poprzednie moduły miały dane parametry planta — m, c, k z modelu MSD; J, R, L z silnika DC; M, m, ℓ z cart-pole. Skąd je wzięliśmy?

  • Z katalogu producenta — najczęściej, ale tolerancje 5–10%.
  • Z pomiarów fizycznych — masa wagą, długość linijką. Dokładne dla mechanicznych, niemożliwe dla elektrycznych „wewnętrznych" (R uzwojenia, L pod napięciem) i nieliniowych (tarcie statyczne).
  • Z first principles — wyprowadzenie z praw fizyki. Wymaga pełnej znajomości geometrii i materiałów.
  • Z identyfikacji — eksperyment + dopasowanie modelu do danych. Standard dla profesjonalnego projektowania kontrolerów.

Identyfikacja to jedyna metoda, która pozwala uchwycić:

  • Rzeczywiste tarcia, lufy, niedokładności — to czego idealne wzory nie modelują.
  • Zmiany w czasie — np. opór zwiększa się ze temperaturą, łożyska zużywają się.
  • Niemodelowane dynamiki — np. rezonanse konstrukcji, których nie ma w idealnym modelu drugiego rzędu.

Z perspektywy kursu: identyfikacja to pierwszy krok projektu sterowania. Bez modelu nie ma LQR/MPC/SMC.

2. Odpowiedź skokowa + fitowanie

Analogia. Wyobraź sobie, że masz tłumik hydrauliczny i chcesz wyznaczyć jego sztywność. Pociągasz za drążek skokowo (zadajesz step), patrzysz jak urządzenie wraca do równowagi (mierzysz y(t)). Z kształtu krzywej (jak szybko zbiega, czy oscyluje) odczytujesz parametry. Komputer robi to systematycznie przez dopasowanie matematycznego modelu.

Najprostszy eksperyment identyfikacyjny:

  1. Wzbudź plant skokiem (z 0 do U_step w t = 0).
  2. Mierz y(t) przez czas dostateczny do osiągnięcia stanu ustalonego.
  3. Dopasuj parametry teoretycznego modelu tak, by jego odpowiedź skokowa minimalizowała SSE względem pomiaru.

Dla planta drugiego rzędu typu MSD parametry to (m, c, k). W demie poniżej możesz spróbować ręcznie — automatyczny optimizer zrobiłby to samo, ale szybciej.

parametry modelu (zgadnij!)

m_fit (masa)1.00
c_fit (tłumienie)2.00
k_fit (sztywność)4.0
σ pomiaru0.040
jakość dopasowania
SSE0.4683R²0.9604
-0.30.180.651.131.601.22.43.64.86pomiary (z szumem)y_clean (true, ukryte!)twój modelt [s]y

Identyfikacja parametryczna to proces dobierania parametrów modelu do zmierzonych danych. W tej wersji wizualnej — ręcznie, suwakami. W praktyce — automatycznie przez optymalizację (najmniejsze kwadraty, RLS, maximum likelihood).

Sum of Squared Errors (SSE):

SSE(θ)=∑i=1N(yimeas−yimodel(θ))2\mathrm{SSE}(\theta) = \sum_{i=1}^N \big(y_i^{meas} - y_i^{model}(\theta)\big)^2

Minimalizacja SSE jako funkcji parametrów θ=(m,c,k)\theta = (m, c, k) daje najlepsze dopasowanie metodą najmniejszych kwadratów. Dla układu liniowego w parametrach (jak nasz MSD) jest to convex problem z analitycznym rozwiązaniem (równania normalne).

R² (coefficient of determination):

R2=1−SSESST,SST=∑(yi−yˉ)2R^2 = 1 - \frac{\mathrm{SSE}}{\mathrm{SST}}, \quad \mathrm{SST} = \sum (y_i - \bar y)^2

Mówi jaką frakcję wariancji w danych wyjaśnia model. R² = 1 = perfekcyjne dopasowanie. R² < 0.5 = model zły lub dane bardzo zaszumione.

Trik: każda kombinacja (m, c, k) różniąca się wzmacniaczem mnożnikowym daje identyczną odpowiedź skokową — można odzyskać tylko stosunki c/m i k/m, nie absolutne wartości. Stąd zwykle ustala się jeden parametr (np. k z prawa Hooke'a) i fitujemy pozostałe.

3. Modele ARX i ARMAX

Modele ciągłe wymagają znajomości struktury fizycznej (jak MSD). W praktyce częściej używa się modeli dyskretnych ARX/ARMAX — czysto matematycznych relacji wejścia-wyjścia, niezależnych od konkretnej fizyki.

ARX (AutoRegressive with eXogenous input)

yk=a1yk−1+a2yk−2+…+anyk−n+b1uk−1+…+bmuk−m+eky_k = a_1 y_{k-1} + a_2 y_{k-2} + \ldots + a_n y_{k-n} + b_1 u_{k-1} + \ldots + b_m u_{k-m} + e_k

Wyjście jest liniową kombinacją swoich poprzednich wartości i wejść. Parametry θ=[a1,…,an,b1,…,bm]\theta = [a_1, \ldots, a_n, b_1, \ldots, b_m] wyznacza się metodą najmniejszych kwadratów:

θ^=(ΦTΦ)−1ΦTY\hat\theta = (\Phi^T \Phi)^{-1} \Phi^T Y

gdzie Φ\Phi to macierz regresji (każdy wiersz to φk=[yk−1,…,yk−n,uk−1,…,uk−m]T\varphi_k = [y_{k-1}, \ldots, y_{k-n}, u_{k-1}, \ldots, u_{k-m}]^T), a YY wektor pomiarów. Klasyczny problem algebry liniowej — dla dobrze postawionej macierzy ΦTΦ\Phi^T\Phi jest jednoznaczne rozwiązanie.

ARMAX (z modelem szumu)

ARX zakłada, że eke_k to biały szum (niezależny w czasie). W praktyce szum często ma strukturę kolorową (skorelowany). ARMAX dodaje model szumu jako średnią ruchomą (MA):

yk=∑aiyk−i+∑bjuk−j+ek+c1ek−1+…+cpek−py_k = \sum a_i y_{k-i} + \sum b_j u_{k-j} + e_k + c_1 e_{k-1} + \ldots + c_p e_{k-p}

Estymacja ARMAX wymaga iteracyjnej procedury (LS nie wystarcza, bo regresory eke_k nie są obserwowalne bezpośrednio). Klasyczne algorytmy: predict error minimization, instrumental variables.

OE, BJ — modele alternatywne

  • OE (Output Error) — y=G(z)u+ey = G(z)u + e bez modelu szumu w równaniu y. Lepsze gdy szum jest tylko pomiarowy, nie procesowy.
  • BJ (Box-Jenkins) — najbardziej ogólny: y=G(z)u+H(z)ey = G(z)u + H(z)e z osobnymi transmitancjami dla planta i szumu.

Wybór modelu zależy od fizyki problemu i jest części sztuki identyfikacji.

4. Recursive Least Squares — online

Analogia. Klasyczne LS to „ważenie z taśmy produkcyjnej": zbierasz wszystkie ważone produkty, dopiero po fakcie obliczasz średnią. RLS to „ważenie ciągłe": z każdym kolejnym produktem aktualizujesz średnią, wykorzystując tylko poprzedni wynik + nowy pomiar — bez powtarzania całej historii. Krytyczne dla embedded i real-time identyfikacji.

Recursive Least Squares przekształca offline algorytm LS w iteracyjną formę:

θ^k=θ^k−1+Kk (yk−φkTθ^k−1)\hat\theta_k = \hat\theta_{k-1} + K_k\,(y_k - \varphi_k^T \hat\theta_{k-1})

Każdy nowy pomiar aktualizuje estymatę o poprawkę proporcjonalną do błędu predykcji (yk−φkTθ^k−1)(y_k - \varphi_k^T \hat\theta_{k-1}). Wzmocnienie KkK_k dynamicznie się dostosowuje na podstawie macierzy kowariancji PkP_k:

Kk=Pk−1 φkλ+φkTPk−1 φk,Pk=1λ(Pk−1−Kk φkTPk−1)K_k = \frac{P_{k-1}\,\varphi_k}{\lambda + \varphi_k^T P_{k-1}\,\varphi_k}, \quad P_k = \frac{1}{\lambda}\left(P_{k-1} - K_k\,\varphi_k^T P_{k-1}\right)

Plus offline LS: nie wymaga pamiętania całej historii pomiarów; konstant memory O(n²) dla n parametrów, constant time per sample. Idealne dla embedded.

Forgetting factor λ∈(0,1]\lambda \in (0, 1] umożliwia śledzenie parametrów zmieniających się w czasie. Dla λ<1\lambda < 1 stare dane są „zapomniane" eksponencjalnie. Trade-off: niższa wartość → szybsze śledzenie, ale wyższa wariancja.

sygnał wzbudzenia

parametry RLS

σ szumu0.05
λ (forgetting factor)0.990
estymaty końcowe (t = 30s)
a_true0.70â0.690b_true0.30b̂0.310
-0.5-0.080.350.771.20612182430â (estymata)a_trueb̂ (estymata)b_truet [s]estymaty parametrów-1.5-0.370.751.8830612182430u (sygnał wzbudzenia)y (zmierzone, z szumem)t [s]sygnały

Recursive Least Squares (RLS) — online wersja najmniejszych kwadratów. Każdy nowy pomiar aktualizuje estymatę parametrów bez ponownego całkowania wszystkich danych.

θ^k=θ^k−1+Kk (yk−φkTθ^k−1)\hat\theta_k = \hat\theta_{k-1} + K_k\,(y_k - \varphi_k^T \hat\theta_{k-1})
Kk=Pk−1φkλ+φkTPk−1φk,Pk=1λ(Pk−1−KkφkTPk−1)K_k = \frac{P_{k-1}\varphi_k}{\lambda + \varphi_k^T P_{k-1}\varphi_k}, \quad P_k = \frac{1}{\lambda}\Big(P_{k-1} - K_k \varphi_k^T P_{k-1}\Big)

Gdzie φk=[yk−1,uk−1]T\varphi_k = [y_{k-1}, u_{k-1}]^T to wektor regresji, θ^=[a^,b^]T\hat\theta = [\hat a, \hat b]^T — estymaty parametrów, λ\lambda — forgetting factor (typowo 0.95–0.99) — daje większą wagę nowym danym.

Persistent Excitation

RLS konverguje do prawdziwych parametrów tylko jeśli sygnał wzbudzenia jest persistently exciting — zawiera „dość dużo" niezależnych komponentów. Dla 2 parametrów potrzebne co najmniej 2 niezależne częstotliwości.

  • Step — tylko jedna informacja (DC + transient). Estymata nie zbiega do prawdy: â i b̂ mogą się stabilizować na innych wartościach, które dają identyczne y dla tego konkretnego wzbudzenia.
  • PRBS — pseudo-random binary. Szerokie pasmo częstotliwościowe, ideał dla identyfikacji. Estymaty zbiegają dokładnie do prawdy.
  • Multi-sin — kilka sinusoid o różnych ω. Też zbiega, ale tylko częstotliwości w sygnale są „widziane".

Spróbuj: tryb „step" — â stabilizuje się gdzieś między 0.5 a 0.9, ale niekoniecznie na 0.7. Przełącz na PRBS — obie estymaty szybko zbiegają do prawdy (0.7, 0.3).

Forgetting factor λ

λ=1\lambda = 1: pamiętaj wszystkie dane (klasyczne LS). Dobre dla parametrów stałych, słabe dla zmieniających się w czasie. λ<1\lambda < 1: starsze dane „zapomniane" — adaptive tracking parametrów zmiennych. Cena: wyższa wariancja estymaty (filtr „zapomnial" historyczne, więc opiera się na mniej danych).

5. Persistent Excitation — warunek istnienia rozwiązania

Analogia. Próbujesz wyznaczyć opór elektryczny R używając prawa Ohma U = R·I. Jeśli mierzysz tylko przy jednym natężeniu prądu (np. I = 1A), masz tylko jedno równanie z jednej niewiadomej — OK. Ale dla R i L (oporu i indukcyjności) potrzebujesz dwóch niezależnych eksperymentów: przy stałym prądzie (mierzy R) i przy zmiennym (mierzy L). Persistent excitation to formalna wersja: ile niezależnych „testów" Twoje sygnały dostarczają.

Definicja

Sygnał uku_k jest persistently exciting rzędu n jeśli istnieje stała α>0\alpha > 0 i N takie że:

1N∑k=1NφkφkT⪰αIndla wystarczająco duz˙ego N.\frac{1}{N}\sum_{k=1}^N \varphi_k \varphi_k^T \succeq \alpha I_n \quad \text{dla wystarczająco dużego } N.

Innymi słowy: średnia z iloczynów zewnętrznych regresorów φk\varphi_k jest dodatnio określona — sygnały „pokrywają" cały podprzestrzeń Rn\mathbb{R}^n.

Reguła praktyczna

Dla n parametrów do estymacji, potrzeba sygnału z n niezależnymi częstotliwościami.

  • Step — n = 1 (DC), zwykle za mało.
  • Sin(ω₁t) — n = 1 (jedna częstotliwość).
  • Sin(ω₁t) + Sin(ω₂t) — n = 2.
  • PRBS — n efektywne wielkie (szerokie pasmo).
  • Chirp — sweep częstotliwości, n „nieskończone" w pewnym paśmie.

Konsekwencje braku PE

Gdy PE nie zachodzi, macierz informacji (ΦTΦ)(\Phi^T\Phi) jest singularna (rank deficient). Klasyczne LS nie ma jednoznacznego rozwiązania — istnieje cała hyperpłaszczyzna parametrów dających ten sam SSE. Wybór konkretnego punktu zależy od stanu początkowego algorytmu albo od regularyzacji.

W RLS bez PE estymata nie zbiega do prawdy — może oscylować lub stabilizować się na innym minimum SSE. Pokazaliśmy to w dema §4 — tryb „step" daje stabilną estymatę, która jest błędna.

6. Identyfikacja w pętli zamkniętej

Wszystko powyżej zakładało, że eksperyment przeprowadzasz przed uruchomieniem regulatora — z plantem w otwartej pętli. W praktyce często to niemożliwe:

  • Plant niestabilny (cart-pole) — otwarta pętla = upadek, żadnych pomiarów.
  • Plant z ograniczonym zakresem operacyjnym — trzymanie w otwartej pętli go niszczy (np. silnik bez sprzężenia pozycyjnego).
  • Plant w produkcji — nie można wyłączyć regulatora dla identyfikacji.

Closed-loop identification radzi sobie z tym, ale ma swoje pułapki:

Problem korelacji

W pętli zamkniętej uk=−K yku_k = -K\,y_k, więc uu jest funkcją yy. Konsekwencja: zaburzenie procesowe ww wpływa na yy, który wpływa na uu, który dalej wpływa na yy. Klasyczne LS „widzi" tę korelację jako fałszywą dynamikę i daje biased estymaty.

Direct vs Indirect Approach

  • Direct: ignoruje regulator, robi LS z (u,y)(u, y) jakby były niezależne. Pracuje gdy szum pomiarowy jest dominujący nad procesowym. Najprostszy w praktyce.
  • Indirect: estymuje cały układ zamknięty (transmitancję od referencji do wyjścia), potem rekonstruuje plant przez znany regulator. Wymaga znajomości regulatora.
  • Joint Input-Output: pracuje z parą sygnałów jako wektorem stochastycznym. Najbardziej rygorystyczne.

External excitation

Najczęściej praktykowane: dodaj zewnętrzny sygnał wzbudzenia r_k (mała sinusoida lub PRBS) do referencji r=rsetpoint+rexcr = r_{setpoint} + r_{exc}. Sygnał wzbudzenia musi być niezależny statystycznie od zaburzeń procesowych — dzięki temu daje informację o plant, nie o regulatorze. Amplituda rexcr_{exc} powinna być dostatecznie mała żeby nie zaburzała operacji, ale duża żeby dawała measurable signal-to-noise.

7. Praktyczny workflow identyfikacji

Profesjonalny workflow identyfikacji parametrycznej:

  1. Analiza fizyki — wybierz strukturę modelu (rząd, liniowość, czy jest opóźnienie). To wymaga inżynierskiej intuicji o procesie.
  2. Projekt eksperymentu — wybierz sygnał wzbudzenia (step, PRBS, chirp), amplitude, częstotliwość próbkowania. Często stosuje się PRBS o pasmie 10× pasmo planta.
  3. Pomiar — uruchom eksperyment, zbierz dane (uk,yk)(u_k, y_k) przez czas wystarczający (typowo kilka razy dłuższy niż dominująca stała czasowa, dla PE = wiele okresów wzbudzenia).
  4. Wstępna analiza — sprawdź jakość danych: czy nie ma saturacji aktuatora, czy szum jest rozsądny, czy nie ma outlierów.
  5. Estymacja parametrów — uruchom LS, RLS, lub iteracyjny estymator (PEM dla ARMAX) i odczytaj θ^\hat\theta.
  6. Walidacja — porównaj z innym zestawem danych (innym wzbudzeniem) niż użyty do estymacji. Sprawdź residual analysis: czy e^k\hat e_k wygląda jak biały szum (autocorrelation, widmo). Jeśli ma strukturę, model jest za prosty.
  7. Iteracja — jeśli walidacja słaba, spróbuj wyższego rzędu modelu, innego typu (ARMAX zamiast ARX), innego sygnału wzbudzenia.

Pułapki praktyczne

  • Brak wzbudzenia — eksperyment „spokojny" nie identyfikuje. Trzeba aktywnie wzbudzać.
  • Saturacja — gdy amplituda wzbudzenia zbyt duża, aktuator wpada w nasycenie. Dane stają się nieliniowe, LS daje biased estymaty.
  • Zmienność w czasie — gdy plant się zmienia w trakcie eksperymentu (np. ogrzewanie się silnika), długie eksperymenty są problematyczne. Lepiej RLS z forgetting.
  • Niemodelowane nieliniowości — fitujesz liniowy model do nieliniowego planta. Estymata jest lokalna wokół punktu pracy. Dla innych warunków może być bardzo nietrafna.
  • Overfitting — zbyt wysoki rząd modelu pasuje doskonale do szumu, ale generalizuje źle. Standardowe lekarstwa: AIC, BIC, cross-validation.

Kanon: Ljung „System Identification: Theory for the User" (1999) — biblia identyfikacji. Söderström & Stoica „System Identification" — bardziej teoretyczne.

Co świadomie pomijam

  • Identyfikacja nieliniowa — Wiener, Hammerstein models, neural networks. Większa złożoność, wykraczająca poza ten moduł.
  • Subspace methods (N4SID, MOESP) — identyfikacja state-space MIMO bezpośrednio bez parametryzacji ARMAX.
  • Frequency-domain identification — używanie odpowiedzi częstotliwościowej zamiast czasowej. Lepsza dla długiego pasma.
  • Maximum Likelihood — bardziej rygorystyczna estymacja przy znanym rozkładzie szumu.