Komputerowa identyfikacja obiektów
ZAJ
CIA V
Zasada najmniejszej sumy kwadratów
" Skąd taki pomysł na dopasowanie ?
" Wyprowadzenie skalarne
" Wersja macierzowa
" Modele nieliniowe względem wejść
" Uogólnienie  aproksymacja z wybraną bazą
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
WPROWADZENIE
Dopasowywanie odpowiedzi modelu do pomiarów poprzez taki dobór parametrów, który daje najmniejszą odległość
odpowiedzi wydaje się naturalnym sposobem wyznaczania parametrów modelu dla obiektu fizycznego. Ponieważ
zmierzone wyjście obiektu jest zmodyfikowanym statycznie lub dynamicznie wejściem obiektu, gdzie sposób
modyfikacji jest opisany z użyciem parametrów, to identyfikacja może być widziana jako działanie odwrotne. Tzn.
na podstawie sygnałów po modyfikacji poszukujemy parametrów tej modyfikacji. Jako miary odległości pomiędzy
sygnałami najczęściej używa się metryki euklidesowej, czyli sumy kwadratów odległości odpowiedzi dla
poszczególnych pomiarów.
Zasada dopasowywania odpowiedzi modelu do zmierzonej w N punktach odpowiedzi obiektu według zasady
najmniejszej sumy kwadratów różnic tych odpowiedzi:
N
2
J1 = w
i - yi
()
"
i =1
nie jest jedynym możliwym wyborem. W pierwszych latach XIX-go stulecia, w czasach gdy Gauss obliczał orbity
planet dopasowując pozycje modelowe do zmierzonych (astronomia była wtedy motorem nauk technicznych),
toczył się spór o to jaka zasada jest właściwa, czy stosować metrykę modułową czy kwadratową. Dzisiaj już wiemy,
że wybór zasady dopasowania zależy od przyjętego kryterium dobrego dopasowania i od rozkładu błędu
pomiarowego. Jednak dziewiętnastowieczna popularność zasady najmniejszej sumy kwadratów okazała się
zgodna z dzisiejszymi, popartymi teorią największej wiarygodności, wyborami dla najczęstszego przypadku
zakłóceń o rozkładzie normalnym.
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
WPROWADZENIE C.D.
Poza teoretycznym uzasadnieniem, estymator najmniejszej sumy kwadratów (ang. Least Squares, w skrócie LS)
ma jedną istotną cechę decydującą o jego powszechnym użyciu  jest wyrażony w formie zamkniętej, czyli jest
dany w postaci wzoru. Wyznaczenie wartości estymat następuje na drodze obliczeń algebraicznych, bez
konieczności stosowania kosztownego czasowo postępowania iteracyjnego. Nie jest tak w przypadku estymatorów
opartych o zasadę najmniejszej sumy modułów:
N
J2 = w
i - yi ,
"
i =1
czy zasadę najmniejszej maksymalnej odległości odpowiedzi
J3 = max w
i - yi .
( )
W przypadku kryterium jakości J1 ustalenie wartości parametrów minimalizujących to kryterium jest możliwe
poprzez obliczenie pochodnych cząstkowych rzędu pierwszego i drugiego względem parametrów. Określenie
analitycznej zależności na minimum kryteriów jakości J2 oraz J3 tą drogą nie jest możliwe.
W dalszej części przypomniano wyprowadzenie postaci estymatora LS dla przypadku obliczeń skalarnych i
przypadku obliczeń macierzowych. Pokazano również zastosowanie estymatora LS w szczególnych przypadkach.
Uwaga:
W dalszej części dzisiejszych zajęć będziemy się zajmować liniowym estymatorem najmniejszej sumy kwadratów,
tzn. takim, który wyznacza parametry modelu liniowego względem wyznaczanych parametrów. Na póz
niejszych
zajęciach zajmiemy się przypadkiem modelu nieliniowego względem parametrów.
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
ESTYMATOR NAJMNIEJSZEJ SUMY KWADRATÓW - RACHUNEK SKALARNY
Dla podanych (zmierzonych) par wartości ui,yi , i = 1,& ,N poszukujemy modelu opisującego zależność pomiędzy
{ }
nimi przy i-tym pomiarze w postaci yi = aui + b + �i , gdzie � jest zakłóceniem pomiarowym, u wielkością wejściową,
Ć
a i b parametrami modelu. Aby otrzymać oszacowania �b parametrów modelu, które w jak najmniejszym stopniu
,
różnią się od rzeczywistych wartości parametrów, można w procesie estymacji zastosować zasadę najmniejszej
Ć
sumy kwadratów odległości pomiędzy wyjściem zmierzonym yi i wyjściem modelu w
i = �ui + b .
NN
2
2
J1 = w
i - yi = yi - au� - b = J1 �,b
()
() ( )
""ĆĆ
i =1 i =1
Wartości estymowanych parametrów należy dobrać w taki sposób, aby wartość kryterium J1 była minimalna. W tym
celu wyznaczmy minimum kryterium J1 w funkcji parametrów. Ekstremum funkcji dwóch zmiennych znajdujemy
"J1 "J1
obliczając pochodne cząstkowe i przyrównując je do zera ( = 0, = 0 ).
Ć
"�
"b
NN N N
"J1 �#
2
= -2 yi - au� - b = 2ś# - yi + b + a =
()
"u "� "� "u ś# 0
i i i �
ź#
"�
i =1 �# #
i =1 i =1 i =1
NN N
"J1 �#
= -2 yi - au� - b = 2ś# - + Nb + � =
()
""w
"u ś# 0
i �
ź#
Ć
"b
i =1 �# #
i =1 i =1
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
ESTYMATOR NAJMNIEJSZEJ SUMY KWADRATÓW - RACHUNEK SKALARNY C.D.
Ć
Następnie możemy wyznaczyć estymowane parametry � oraz b przekształcając powyższe zależności do postaci:
NN
ż#
Ć
"y = Nb + �"u
ii
�#
�#
i =1 i =1
�#
NN N
2
�#
Ć
"uy = b"u + �"u
i i i i
�#
�# i =1 i =1 i =1
Ć NN N N
ż#
y = �u + b
1
�#
2
u =
�#
"u , y = 1 "y , uy = 1 "u yi, u2 = 1 "u
ii i i
Ć NN N N
i =1 i =1i =1 i =1
�#
�#uy = bu + �u2
Rozwiązując powyższy układ równań ze względu na estymowane parametry otrzymujemy:
Ć
b = y - �u
2
uy = uy - a� + a�2
uy - uy
� =
2
u2 - u
Aby mieć pewność, że parametry przez nas wyznaczone odnoszą się na pewno do minimum kryterium J1 należy
policzyć drugie pochodne tej funkcji.
"2J1 N 2 "2J1 "2J1 "2J1 N
= 2 = 2N == 2
"u "u
�
Ć Ć
"�2 i =1 i
"b2 "�"b "b"�
i =1
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
ESTYMATOR NAJMNIEJSZEJ SUMY KWADRATÓW - RACHUNEK SKALARNY C.D.
Aby odszukane ekstremum było minimum muszą być spełnione dwa następujące warunki:
"2J1 "2J1
2
Ć
�# ś#
"2J1 "2J1 "2J1
"�2
"�"b
1). W => 0 �! W = - > 0
"�2 Ć ś# Ć ź#
"b2 �# "�"b
"2J1 "2J1
#
Ć
Ć
"b"� "b2
NN N
2
W = 4N - 4
"u "u "u > 0
ii i
i -1 i =1 i =1
"2J1 "2J1 N 2
2). > 0 = 2 > 0
"u
"�2 "�2 i =1 i
W analizowanym przez nas przypadku obydwa te warunki są spełnione.
Ćwiczenie:
Pobór mocy z sieci energetycznej w małym mieście ma pewną część niezmienną i część zależną od temperatury
otoczenia. Jako zmienną wejściową przyjmijmy odchyłkę temperatury "T od 25�C. Mierząc pobór mocy P [MW] w
funkcji odchyłki temperatury "T uzyskano następujące wyniki:
0 -2 -5 -10 -15 -20 -25 -30 -35 -40 -45
"T
P 20 22 25 32 40 53 60 69 80 84 90
Czy liniowy model zjawiska jest uzasadniony ? Jaką postać będzie miał model zjawiska ? Jak wyznaczyć wartości
parametrów skalarnym estymatorem LS ?
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
ESTYMATOR NAJMNIEJSZEJ SUMY KWADRATÓW - RACHUNEK MACIERZOWY
Przedmiotem naszego badania jest obiekt wielowymiarowy o n zadanych wejściach u(1),u(2),u(3),...,u(n) i jednym
mierzonym wyjściu y. Pomiarów tego obiektu dokonujemy N-krotnie (i = 1,2,3,..,N ), a wyniki takich pomiarów
zawiera następująca tabela:
( ( ( (
u11) u12) u13) u1n) y1
( ( ( (
u21) u22) u13) u2n) y2

1 3
( ( ( (
uN ) uN2) uN ) uNn) yN
Przyjmujemy, że badany obiekt jest liniowy i pomiary są opisane modelem
yi = a0 + a1ui(1) + a2ui(2) +& + anui(n) + �i gdzie i = 1,2,& ,N
gdzie � jest przypadkowym zakłóceniem pomiarowym.
Oszacowanie (czyli estymatę) �0,�1,& ,�n parametrów a0,a1,& ,an tego modelu znajdziemy minimalizując sumę
kwadratów odchyłek
NN'" 2 N
2
2
J1 = = yi - �0 - �1ui(1) -& - �nui(n) = 4
1 �0,a1,& ,�n
()
()
"e = "�# yi - yi ś# "
ś#
i ź#
�# #
i =1 i =1 i =1
względem parametrów �0,�1,& ,�n .
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
ESTYMATOR NAJMNIEJSZEJ SUMY KWADRATÓW - RACHUNEK MACIERZOWY C.D.
W tym celu obliczamy pochodne cząstkowe funkcji J1 względem tych parametrów, a ich wartości optymalne
wyznaczamy z warunku optymalności, czyli przyrównania pochodnych cząstkowych do zera. Aby jednak
wyznaczenie tych parametrów było łatwiejsze do obliczeń posłużymy się rachunkiem macierzowym. W tym właśnie
celu wprowadzamy następujące oznaczenia:
( ( (
Ą#1 u11) u12)
& u1n) ń#
ó#Ą#
( ( (
1 u21) u22) & u2n)
ó#Ą#
ó#Ą#
( ( (
MACIERZ WEJŚĆ U =
& u3n) Ą# ,
ó#1 u31) u32)
ó# Ą#
ó#Ą#
1
( ( (
ó#1 uN ) uN2) & uNn) Ą#
Ł#Ś#
w której dodano jedynki w pierwszej kolumnie aby ujednolicić zapis macierzowy dla wyrazu wolnego a0 w modelu.
y1 w
1
Ą# ń# Ą# ń#
ó# Ą# ó# Ą#
y2 w
2
ó# Ą# ó# Ą#
WEKTORY WYJŚĆ obiektu i modelu Y = ó# y3 Ą# v
= ó# w
3 Ą#
ó# Ą# ó# Ą#

ó# Ą# ó# Ą#
ó# Ą# ó# Ą#
Ł#yN Ś# Ł#w
N Ś#
e1 �0
Ą# ń# Ą# ń#
ó#e Ą# ó#� Ą#
2 1
ó# Ą# ó# Ą#
Ć
WEKTOR ODCHYA
EK WYJŚĆ obiektu i modelu e = ó#e3 Ą# = v
- Y , WEKTOR ESTYMAT PARAMETRÓW modelu � = ó#�2 Ą#
ó# Ą# ó# Ą#

ó# Ą# ó# Ą#
ó# Ą# ó# Ą#
Ł#eN Ś# Ł#�n Ś#
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
ESTYMATOR NAJMNIEJSZEJ SUMY KWADRATÓW - RACHUNEK MACIERZOWY C.D.
Tak więc modelowa zależność opisująca pomiary przybierze postać macierzową:
Y = U� + �
Minimalizowana suma kwadratów błędów przyjmuje w zapisie macierzowym postać
T
T
Ć
J1 = eTe = Y - U� v
- U� = v
T Y - YTU� - �T�T v
+ �T�U� =
( ) ( )
T
= YT Y - 2�T�T Y + �T�U� = 4
1 �
(Ć)
Ć
Różniczkując teraz funkcję J1 � względem wektora � i przyrównując do zera otrzymujemy warunek optymalności
(Ć)
"J1
= -2UT Y + 2UT۸ = 0
Ć
"�
skąd
T
Ć
UU� = UT Y
Ostatnie równanie jest nazywane równaniem normalnym.
-1
T
Po lewostronnym pomnożeniu obu stron równania normalnego przez UU otrzymujemy
( )
-1
T
Ć
� = UU UT Y
( )
co jest klasycznym wzorem na estymator LS.
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
ESTYMATOR NAJMNIEJSZEJ SUMY KWADRATÓW - RACHUNEK MACIERZOWY C.D.
T T
Aby jednak otrzymać rozwiązanie, macierz UU nie może być osobliwa (det UU `" 0 ). W praktyce oznacza to, iż
( )
liczba pomiarów N musi być co najmniej równa n + 1, czyli ilości parametrów estymowanych, oraz że w macierzy
wejść elementy jednej kolumny nie mogą być liniowo zależne od elementów innych kolumn
"J1
Można też dowieść, że przyjęty przez nas warunek optymalności ( = 0) i wynikająca stąd zależność na
Ć
"�
Ć
wyznaczenie wektora współczynników � prowadzą do minimalizacji a nie maksymalizacji funkcji J1 � , gdyż
(Ć)
"J1 "2J1 T
T
różniczkując jeszcze raz otrzymujemy = 2UU, a z dodatniej określoności macierzy UU wynika dodatnia
Ć Ć
"� "�2
określoność macierzy drugich pochodnych cząstkowych, co stanowi warunek na minimum funkcji J1 � .
(Ć)
Ć
Po wyznaczeniu estymaty � wektora parametrów modelu � , możemy obliczyć wektor wyjść modelu odpowiadający
Ć
macierzy pomiarów v
= ۸, oraz estymatę (predykcję) wyjścia obiektu w
i = ui� przy dowolnej innej wartości wejść,
Ą#1
gdzie ui = ui(1) ui(2) & ui(n) ń# . Inaczej mówiąc, jesteśmy w stanie przewidzieć zachowanie się obiektu w innych
Ł# Ś#
warunkach, przy nie obserwowanej dotąd kombinacji wartości wejściowych. Poprawność takiej predykcji jest
ostatecznym potwierdzeniem poprawności modelu.
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
Przykład: Funkcje Matlaba związane z tematem: operator  \ , funkcja polyfit
Poniżej przedstawiono przykładowy skrypt Matlaba estymujący parametry modelu liniowego y = au + b dwiema
metodami.
% Wartości wielkości wejściowej
u=[0.1:0.1:10]';
% Pierwsza seria pomiarów dla zadanych wartości u
y1=pomiar(u);
% Druga seria pomiarów, będzie inna realizacja zakłóceń
y2=pomiar(u);
subplot(2,1,1)
% Obliczenie parametrów z pierwszej serii pomiarów 15
U=[u,ones(m,1)];
10
A1=U\y1;
5
plot(u,y1,'rx',u,U*A1,'b-');
0
0 2
subplot(2,1,2)
15
% Obliczenie parametrów z drugiej serii pomiarów
A2=polyfit(u,y2,1);
10
plot(u,y2,'rx',u,A2(1)*u+A2(2),'b-');
5
0
Na powyższych wykresach przedstawione zostały efekty estymacji
0 2
parametrów a i b, z zastosowaniem operatora dzielenia lewostronnego,
Ć
który jest interpretowany w sensie najmniejszej sumy kwadratów dla równań nadokreślonych (� = 1,024 b = 1,4605 ,
Ć
górny wykres), oraz z użyciem funkcji polyfit (� = 0,9835 b = 1,5536 , wykres dolny). Wartości wyznaczone
poszczególnymi metodami różnią się, ponieważ dopasowanie wykonano dla różnych realizacji zakłóceń.
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
UWAGI O MODELACH NIELINIOWYCH WZGL
DEM ZMIENNYCH NIEZALEŻNYCH, ALE LINIOWYCH
WZGL
DEM PARAMETRÓW
Zauważmy, że w powyższych wyprowadzeniach określenie model liniowy dotyczy liniowości względem
parametrów, a nie liniowości względem zmiennych niezależnych. I tak współczynniki ai modelu o ogólnej postaci
n
i
y = �" f t , gdzie f(i)(t) są nieliniowymi funkcjami parametru t, mogą być wyznaczone wyprowadzonymi
( )
"a ( )
i
i =1
powyżej estymatorami. W zapisie macierzowym widać podobieństwo do klasycznego sformułowania problemu LS.
1 n
( ) ( )
y1 Ą# f t1 f t1
Ą# ń# ( ) ( )ń#
Ą#
ó# Ą#ó#
= �" �
Ą#
ó# Ą#ó#
n
( )
ó#yN Ą#ó#f 1 tN f tN Ś#
Ł# Ś#Ł# ( )( ) ( )Ą#
Przykład
Problem wyznaczenia współczynników modelu siły oporu ruchu o kwadratowej zależności od prędkości można
rozpatrzyć jako liniowy, przyjmując jako wielkość wejściową kwadrat prędkości v2.
v2=x F
v
F
a �"v2 + b a �" x + b
Tak określony model możemy opisać układem równań :
F1 v12 1
Ą# ń#Ą#ń#
a
ó# Ą#ó# Ą# ń#
= Ą# �"
ó#bĄ#
ó# Ą#ó#Ą#
ó#FN Ą#ó#vN 2 1Ą# Ł# Ś#
Ł# Ś#Ł#Ś#
i estymować parametry a, b estymatorem LS.
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
Przykład: Korekcja zniekształceń kamery cyfrowej
Przeanalizujmy program korygujący zniekształcenia beczkowe zdjęć z aparatu zrobionych w trybie makro.
Opieramy się na pokazanym obok testowym obrazie kratowym. Problem rozwiniemy na ćwiczeniach.
close all
clear all
i1=imread('Krata1.jpg');
i1g=double(rgb2gray(i1));
[sx, sy]=size(i1g);
cx=sx/2;
cy=sy/2;
% wyszukanie przecięć linii
wzorzec=[
0 0 0 1 1 1 0 0 0
0 0 0 1 1 1 0 0 0
0 0 0 1 1 1 0 0 0
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
0 0 0 1 1 1 0 0 0
0 0 0 1 1 1 0 0 0
0 0 0 1 1 1 0 0 0
];
ic=normxcorr2(wzorzec,i1g);
ig=ic<-0.7;
ig=bwmorph(ig,'shrink',Inf);
[ix,iy]=find(ig>0);
r=sqrt((ix-cx).^2+(iy-cy).^2);
[rs,i]=sort(r);
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
% ustalenie punktów odniesienia (siatki regularnej)
x0=ix(i(1)); y0=iy(i(1));
300
dx=abs(ix(i(1))-ix(i(4)))-1;
dy=abs(iy(i(1))-iy(i(4)))-3;
250
jx=x0+dx*[-4:1:5];
jy=y0+dy*[-6:1:7];
200
[jX, jY]=meshgrid(jx,jy);
jx=jX'; jx=jx(:);
150
jy=jY'; jy=jy(:);
rm=sqrt((jx-cx).^2+(jy-cy).^2);
100
[rms,j]=sort(rm);
50
imshow(ig), hold on
plot(cy,cx,'r+') % punkt centralny
0
0 50 100 150 200 250 300
plot(jY,jX,'bo') % punkty odniesienia
% identyfikacja modelu algorytmem LS
% uzupełnij wyliczanie modelu stopnia 3 zależności rs od rms
% wartości wyliczone z modelu zapisz do zmiennej R
???
???
???
figure, plot(rms,rs,'b.', rms,R,'g-', [0 3e2],[0 3e2],'r-')
% prostowanie obrazka
xx=1:sx;
yy=1:sy;
[XX,YY]=meshgrid(xx,yy);
RR=sqrt((XX-cx).^2+(YY-cy).^2);
RRr=RR*A(1)+RR.^2*A(2)+RR.^3*A(3);
RRu=RRr./RR;
i1r=interp2(i1g,(YY-cy).*RRu+cy,(XX-cx).*RRu+cx);
figure, subplot(2,1,1),imagesc(i1g), subplot(2,1,2),imagesc(i1r), colormap gray
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
ZADANIA  LINIOWA ESTYMACJA LS
Zadanie 1
Na belce z mostkiem tensometrycznym (patrz rysunki poniżej) dokonano pomiarów w celu wzorcowania układu do
pomiaru siły. Na zmierzony stosunek napięcia wyjściowego z mostka do napięcia zasilającego ma wpływ
przyłożona do belki siła gnąca F, odchyłka "T temperatury otoczenia T od temperatury znamionowej T0, błąd
zerowania i błąd przypadkowy pomiaru.
R1
R2
F
T
Uzas
R3 R4
U
Możemy więc zdefiniować model układu w postaci następującej:
U
ur = = kFF + kT "T + u0 + �
Uzas
Przeprowadzone pomiary dały wyniki:
F [kN] 1 2 3 4 5 6 7 8 9 10
"T [K] 0.5 2 5 1 1.2 0.2 -1.2 -3 -2 -1
ur [V/V] 0.0124 0.0226 0.0336 0.0419 0.0509 0.0628 0.0705 0.0813 0.0903 0.1014
Należy dokonać estymacji parametrów powyższego modelu macierzową wersją estymatora LS.
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
Zadanie 2
Rozkład sygnału okresowego na szereg Fouriera może być widziany jako problem estymacji parametrów modelu
liniowego względem parametrów, ale nieliniowego względem zmiennej niezależnej, rozwiązywany wg zasady LS.
W uproszczonym przypadku sygnału nieparzystego wyznaczane muszą być tylko amplitudy składowych
sinusoidalnych. Na poniższym schemacie przedstawiono interpretację wejściowo-wyjściową tego przekształcenia,
gdzie � reprezentuje część sygnału, której widmo leży poza n-tą harmoniczną.
a1
sin(�t)
a
sin(2�
f(t)
an
sin(n�
�
Zajmijmy się przypadkiem dyskretnym, czyli algorytmem DFT na podstawie N próbek funkcji f(t) w jednym okresie:
N -1
i
F n = i sin�# 2Ą n
( ) ( )
"f �#ś#
ś#ź#
N
i =0 #
Należy wykazać, że istotnie algorytm DFT jest równoważny (z dokładnością do współczynnika skalującego)
algorytmowi LS zastosowanemu do powyższego równania, oraz przeprowadzić obliczenia dla 100 próbek
nieparzystego przebiegu prostokątnego rozkładanego na 10 harmonicznych.
Katedra Metrologii AGH Kraków 2006
Komputerowa identyfikacja obiektów
Zadanie 3
Uzupełnij algorytm korekcji zniekształceń geometrycznych obrazu z kamery cyfrowej przedstawiony na wykładzie (i
dostępny na serwerze) o formowanie macierzy i obliczenia estymatora LS. Sprawdz
wizualnie jakość korekcji na
obrazie.
LITERATURA DODATKOWA
Sydenham P.H., Podręcznik Metrologii, WKiA
Warszawa 1988 (rozdział 8 pt. Estymacja parametru)
de Larminat P., Thomas Y., Automatyka  układy liniowe, tom 2 Identyfikacja, WNT Warszawa 1983
Mańczak K., Nahorski Z., Komputerowa identyfikacja obiektów dynamicznych, PWN, Warszawa 1983
Bubnicki Z., Identyfikacja obiektów sterowania, PWN, Warszawa 1974
Eykhoff P., Identyfikacja w układach dynamicznych, PWN, Warszawa1980
Katedra Metrologii AGH Kraków 2006