5 - Identifikace
Michael Šebek Automatické řízení 2016 8-3-16
Identifikace Automatické řízení - Kybernetika a robotika
Aneb jak získat model systému z dat (a validovat ho na jiných datech) • white box (víme vše): ze základních principů (fyz-chem-bio-…) • grey box (víme něco): známe třeba typ modelu, hledáme parametry • black box (nevíme nic): neznáme ani typ modelu, řád, nelinearitu, ... V ARI ukážeme jednoduché experimentální metody, off-line, open-loop • z časové odezvy • z frekvenční odezvy • základy nejmenších čtverců V dalších předmětech • sofistikované stochastické metody typu black box • rekurentní, on-line, v uzavřené smyčce • pro pokročilé zájemce (a diskrétní systémy) L. Ljung: System Identification: Theory for the User (2nd Ed.) Prentice Hall, 1999. ISBN 978-0136566953 • Matlab: System Identification Toolbox Michael Šebek
ARI-05-2015
2
Aproximace ze skokové odezvy - Bump test Automatické řízení - Kybernetika a robotika
Hledáme lineární odchylkový model – tomu přizpůsobíme experiment • pozor na pracovní bod a velikost skoku Postup 1. Experimentálně změříme skokovou odezvu → graf nebo tabulka 2. Na ní naměříme několik vybraných bodů 3. Dosadíme hodnoty do obecně vypočítané odezvy a řešíme rovnice pro neznámé parametry – je to obtížné, tak hledáme zvláštní hodnoty Metody • Klasické, jednoduché metody, off-line, open-loop • Deterministické - fungují, jen když nejsou (jsou malé) šumy, statistické vlastnosti neznáme • Soustava musí být stabilní! • Pokud má systém dopravní zpoždění, odečteme ho předem • Pokud má systém ofset, odečteme ho předem Michael Šebek
ARI-05-2015
TD
3
1. Řád bez nuly Automatické řízení - Kybernetika a robotika
y (∞ )
Aplikujeme u ( s )= u (∞) s a hledáme k G (s) = 1 + Ts
Změříme y (∞) a vypočteme k =y (∞) u (∞) Najdeme 0.63 y (∞) a odměříme T
1. 2. Pro 1.
0, 63 y (∞)
G (s) =
u (∞ ) y (1)
y (1) =
Michael Šebek
y (t ) = (1 − e )ku (∞)
0
y (t ) =
Konečně pro je u (∞ ) t Ti
t T
0
1 Ti s
Odměříme y (1) a vypočteme Ti =
−
T G (s) =
1 + Tp s 1 + kp = Ti s Ti s
u (∞ ) Ti
u (∞ ) Ti
k p u (∞ ) 1
1 ARI-05-2015
4
2. řád bez nul - kmitavý případ Automatické řízení - Kybernetika a robotika
Hledáme
ω G (s) = k 2 s + 2ζωn s + ωn2 2 n
1.
%OS y (∞ )
±2%
Změříme
y (∞), %OS , Ts
2.
Vypočteme
ζ =
Michael Šebek
Ts
− ln ( %OS 100 )
y (∞ ) 4 k = , ωn = , 2 2 T u (∞ ) ζ π + ln ( %OS 100 ) s
ARI-05-2015
5
Jiná metoda pro 2. řád bez nul kmitavý Automatické řízení - Kybernetika a robotika
• Hledáme
y (t )
ωn2 G (s) = k 2 s + 2ζωn s + ωn2 • Aplikujeme u (∞ ) u ( s) = s 1. Změříme y (∞), A1 , A2 , Td
2.
Vypočteme = k
A1
y (∞ )
A2
Td
A1 y (∞ ) = , µ ln= ,ζ u (∞ ) A2
µ
= , ωn
4π 2 + µ 2
2π Td 1 − ζ 2
• Klasické názvosloví:
A1 A2 faktor útlumu, T0 časová konstanta µ tzv. logaritmický dekrement útlumu • Pro zajímavost A1 A 1 2π ln 1 = ζωnTd , ωd = = eζω T , µ = dále platí A2 n − 1 An Td n d
Michael Šebek
ARI-05-2015
6
Zvláštní systém integračního charakteru bez nul Automatické řízení - Kybernetika a robotika
Pro systém s „integračním chováním“ G ( s) =
k s (Ts + 1)
a
u (s) =
1 s
k k Tk Tk = − + s 2 (Ts + 1) s 2 s s + 1 T h(t ) =kt − kT + kTe −t T
1. 2. 3.
τ
nakreslíme asymptotu v nekonečnu odečteme T τ odečteme τ a vypočteme k =
T
T
Obdobně pro složitější 1. 2.
hasymptota (t ) =kt − kT =k ( t − T )
G ( s) =
k s (Ts + 1)
n
k = směrnice asymptoty a platí n −1
Michael Šebek
h(T ) n = e− n k (n − 1)!
h(T )
n h(T ) k ARI-05-2015
1
2
T 3
4
5
0.37 0.27 0.22 0.20 0.18 7
Obecnější systém integračního charakteru Automatické řízení - Kybernetika a robotika
Pro systém s „integračním chováním“ G (s)
1 1 u (∞ ) , ( ) u s = TI s n s 1 T s + ( ) ∏ k
u (∞ )
k =1
dy (t ) u (∞ ) lim s 2 y ( s ) ) lim ( sG (= s )u (∞) ) = lim= ( t →∞ dt s →0 s →0 TI
1. 2. 3. 4.
Nakreslíme asymptotu v nekonečnu Odečteme její směrnici a vypočteme (odečteme) TI Potom uděláme derivaci odezvy a z ní identifikujeme systém proporcionálního charakteru
Téhož dosáhneme použitím impulzního vstupu u (t ) = cδ (t ) u (s) = c
Michael Šebek
yimpulse = ( s ) G= ( s )c GP ( s ) c s ARI-05-2015
TI
( s ) sG ( s ) y= sy = D (s) G= sG (s) = P (s)
1 TI
u (∞ ) s 1
n
∏ (T s + 1) k =1
k
Jak realizovat Dirac? Krátkým obdélníkovým pulsem s plochou c ! 8
2. řád bez nul - nekmitavý případ Automatické řízení - Kybernetika a robotika
Pro G ( s)
y
k e − sL , T2 ≤ T1 (1 + T1s )(1 + T2 s )
T2 T1 ∈ [ 0.1,1]
je skoková odezva T1 T2 − ( t − L ) T1 − ( t − L ) T2 1 e e − − k T1 ≠ T2 T −T T2 − T1 1 2 y (t ) = k 1 − e − (t − L ) T1 − t e − (t − L ) T1 T1 = T2 T 1
y • V horním grafu odezev pro T2 T1 ∈ [ 0.1,1] ≈ 0.72 se zdají být různé, ale po normalizaci času na t (T1 + T2 ) je vidět, že jsou podobné • A proto je těžké určit parametry robustně ze skokové odezvy, šlo by to lépe z impulzní • Protínají se přibližně v bodě(τ ≈ 1.256, f (τ ) ≈ 0.72 ), toho využívá Strejcova metoda – viz příklady Michael Šebek
ARI-05-2015
t
T2 T1 ∈ [ 0.1,1]
τ ≈ 1.256
= τ t (T1 + T2 ) 9
Systém derivačního charakteru (s nulou v nule) Automatické řízení - Kybernetika a robotika
• Pro systém „derivačního charakteru“ (tedy s nulou v s = 0) G(s) =
totéž „po integraci“
TD s ∏ (Ti s + 1)
u (∞ ) = (s) y ( s ) G= s
Skoková odezva derivačního členu
TD u (∞ ) ∏ (Ti s + 1)
y (∞) =0
TD u (∞ ) u (∞ ) = yI ( s ) = y ( s ) s ∏ (Ti s + 1) s
1. 2.
Skokovou odezvu derivačního členu nejprve integrujeme A pak identifikujeme vhodnou z předchozích metod
ct ⇔ u ( s ) = c s2 . • Téhož dosáhneme vybuzením rampou u (t ) = Pak výstup rovnou odpovídá skokové odezvě c ( ) ( ) y= s G = s systému „proporcionálního charakteru“ ramp 2 s
Michael Šebek
ARI-05-2015
TD c ∏ (Ti s + 1) s 10
Identifikace z frekvenční odezvy Automatické řízení - Kybernetika a robotika
• Odezva většinou naměřená (analyzátorem, typicky 0-100kHz), ale • může být i „vypočtená“ (příklad: mechanický model a metoda konečných prvků) Metody • Podívat se, odhadnout vlastnosti a zkusmo napasovat asymptoty (Nise paperback 10.13 strana 665) • Obecné metody interpolace, fitting, nejmenší čtverce • Speciální metody (starší) pro Bodeho nebo Nyquistův graf Michael Šebek
ARI-05-2015
11
Nejmenší čtverce – Least Squares Automatické řízení - Kybernetika a robotika
• Přeurčená soustava lineárních rovnic ( A je m × n, m > n )
Ax = b pokud rank A ≠ rank [ A
b ] nemá řešení!
• Varianta - Nejmenší čtverce
min Ax − b = x
∑ (∑ m
n
=i 1 =j 1
aij x j − bi
)
2
• ekvivalentní minimalizaci bez odmocniny (kvadrátu normy) • = r Ax − b se nazývá reziduum nebo odchylka • Řešení minimalizující normu rezidua se nazývá řešení s nejmenšími čtverci • Pro A plné sloupcové hodnosti najdeme řešení pomocí pseudoinverze:
x = ( A A ) AT b T
Michael Šebek
ARI-05-2016
−1
12
Data fitting Automatické řízení - Kybernetika a robotika
• Vhodným výběrem koeficientů
x1 , x2 , , xn
• napasujte funkci (lineární kombinaci bázových funkcí neboli regresorů gi ( t ) ) g= ( t ) x1 g1 ( t ) + x2 g 2 ( t ) + + xn g n ( t )
• na data (neboli měření) • tak, aby
( t1 , y1 ) , ( t2 , y2 ) , , ( tn , yn )
g ( t1 ) ≈ y1 , g ( t2 ) ≈ y2 , , g ( tn ) ≈ yn
= g ( t1 ) y= y2 ,= , g ( t n ) yn • Obvykle je m n a neexistuje přesné řešení 1 , g ( t2 )
• Takže hledáme nejlepší řešení ve smyslu nejmenších čtverců
min ∑ i =1 ( x1 g1 ( ti ) + x2 g 2 ( ti ) + + xn g n ( ti ) − yi ) m
Michael Šebek
ARI-05-2016
2
13
Data fitting Automatické řízení - Kybernetika a robotika
• Data fitting převedeme na maticový problém
min Ax − b x
• pomocí
g1 ( t1 ) g 2 ( ti ) g1 ( t2 ) g 2 ( t2 ) = A g1 ( tm ) g 2 ( tm )
Michael Šebek
g n ( t1 ) g n ( t2 ) ,x = g n ( tm )
ARI-05-2016
x1 x 2 = ,b xn
y1 y 2 ym
14
Data fitting s polynomy Automatické řízení - Kybernetika a robotika
• Pro
g ( t ) =x1 + x2t + x3l 2 + xnt n −1
• Jsou bázové funkce • a
1 t1 1 t2 A = 1 tm
k −1 g= t t , k − 1, , n ( ) k
t12 t1n −1 x1 x t22 t2n −1 2 = , x = ,b n −1 2 tm tm xn
y1 y 2 ym
• Při interpolaci je m = n a g ( ti ) = yi splníme přesně řešením Ax = b • Při aproximaci je m > n a snažíme se o malou odchylku Michael Šebek
ARI-05-2016
min Ax − b x
15
Identifikace metodou nejmenších čtverců Automatické řízení - Kybernetika a robotika
• Diskrétní lineární model, nazývaný v oblasti identifikace Auto-Regressive model with eXogenous input (ARX) • Daný buď lineární diferenční rovnicí se stochastickým členem y ( t ) + a1 y ( t − 1) + + an y ( t − n= b0u ( t ) + b1u ( t − 1) + + bnb u ( t − nb ) + e ( t ) a)
• nebo přenosem (polynomiálním popisem) a ( d ) y (t ) = b ( d ) u (t ) + e (t )
b(d ) 1 y (t ) = u (t ) + e (t ) a (d ) a (d )
a ( d ) =1 + a1d + + ana d na , b ( d ) =b0 + b1d + + bnb d nb
• Pokud je předem známe pevné dopravní zpoždění a ( d ) y= ( t ) b ( d ) u ( t − Td ) + e ( t )
Michael Šebek
y= (t )
ARI-05-2016
b(d ) 1 u ( t − Td ) + e (t ) a (d ) a (d )
16
Jednorázová identifikace Automatické řízení - Kybernetika a robotika
• Celý soubor naměřených dat se zpracuje najednou • Kompaktně zapsáno = b Ax + r , kde − y ( −1) − y (1 − na ) u (1) u ( 0) y (1) − y ( 0) − y ( 2 − na ) u ( 2 ) y − y ( 0) u (1) 2 ( ) A − y (1) b = y ( m ) − y ( m − 1) − y ( m − 2 ) − y ( m − na ) u ( m ) u ( m − 1) a1 e (1) ana e 2 ( ) x = r b0 e m ( ) bnb
u (1 − nb ) u ( 2 − nb ) u ( m − nb )
je vektor měření výstupů, matice dat, vektor parametrů a vektor chyb predikce, • Hledáme min= r Ax − b x
Michael Šebek
ARI-05-2016
17
LS identifikace – další jemnosti Automatické řízení - Kybernetika a robotika
• Stochastický (Bayesovský) přístup – důkazy • Numerická implementace • jednorázová identifikace • průběžná identifikace - rekurzivní postup • Proměnné parametry • zapomínání, směrové zapomínání • adaptivní řízení • Zabudování apriorní informace • Identifikovaný systém není dostatečně vybuzen • lineární závislost dat – např. identifikace v uzavřené smyčce • návrh experimentu / volba budicího signálu Michael Šebek
ARI-05-2016
To vše v dalších předmětech 18
