Proceedings of International Scientific Conference of FME Session 4: Automation Control and Applied Informatics
Paper 27
Srovnání přístupů zjišťování kauzality, statistické analýzy a identifikace dynamických soustav MORÁVKA, Jan1 1
Ing., Ph.D., TŘINECKÉ ŽELEZÁRNY inženýring, a.s., Středisko elektro, MaR a ASŘ, 739 70 Třinec,
[email protected], http://www.trz.cz
Abstrakt: Statistická analýza a identifikace systémů v průmyslové praxi se provádí hlavně s cílem zjištění účinku předpokládaných akčních veličin pro následný návrh a realizaci řídicích (sub)systémů. Obecněji jde o zjištění kauzálních i kvalitativně-kvantitativních vztahů mezi veličinami technologických procesů a získání matematických modelů soustav. V příspěvku jsou srovnávány různé přístupy k této úloze: expertní (vizuální), klasický statistický, teorie náhodných procesů, analýzy časových řad, moderní statistický a identifikace na základě Box-Jenkinsovy metodologie. Přístupy jsou dokumentovány na příkladu stacionárních veličin (kontinuálních procesů) vstupu a výstupu proporcionální soustavy se setrvačností 1.řádu. Při povrchní analýze přitom nastává jev tzv. zdánlivé nezávislosti (nekorelovanosti, nonkorelace) veličin, který je duálním jevem ke zdánlivé závislosti (korelaci). Klíčová slova: identifikace, statistická analýza, kauzalita, dynamické systémy
1 Úvod Analýza složitých technologických procesů (mezi něž patří např. procesy hutnické, zvláště pak proces vysokopecní a aglomerační) si kromě komplexního popisu klade za cíl: • zjistit kauzální závislosti a kvalitativně-kvantitativní vztahy mezi (měřenými) veličinami • určit vhodné akční a pomocné akční veličiny pro realizaci řídicích (sub)systémů • získat matematický model soustavy (strukturální a parametrická identifikace). Vzhledem k vysoké náročnosti uvedeného problému je nutné používat komplexní analytický přístup, který kombinuje vzájemně se prolínající metody: ! klasické metody matematické statistiky (vícerozměrná regrese průřezových dat, testování regresního tripletu, analýza reziduí) [ANDĚL, J. 1985], [MELOUN, M. & MILITKÝ, J. 1994], [WONNACOT, T. H. & WONNACOT, R. J. 1993], [HEBÁK, P. & HUSTOPECKÝ, J. 1987] ! statistickou dynamiku (teorie náhodných procesů, korelační a spektrální metody) [BENEŠ, J. 1961], [TŮMA, J. 1997] ! statistickou analýzu časových řad (Box-Jenkinsova metodologie, korelační a spektrální analýza) [ANDĚL, J. 1975], [CIPRA, T. 1986] ! moderní metody matematické statistiky a analýzy časových řad (nestacionární a vektorové procesy, kointegrace, Grangerova kauzalita, analýza impuls-reakce, kotrendová analýza atd.) [ARLT, J. 1999], [BIERENS, H. 1999a,b,c], [JOHANSEN, S. 1997], [VÍŠEK, J. Á. 1998] ! ekonometrii (soustavy simultánních dynamických rovnic) [HUŠEK, R. & WALTER, J. 1976], [GARAJ, V. & ŠUJAN, I. 1980] ! modelování, simulaci a identifikaci dynamických systémů [LJUNG, L. 1987], [ZÍTEK, P. 1990], [NOSKIEVIČ, P. 1999], [TŮMA, J. 1998] ! teorii automatického řízení (diskrétní modely spojitých soustav) [VÍTEČEK, A. 1988], [VÍTEČEK, A. & WAWRZICZKOVÁ, M. 1988].
V průmyslové praxi je úloha analýzy a identifikace navíc ztížena nedostatečnosti teoretického popisu zkoumaného technologického procesu, absencí a problematičnosti měření některých relevantních veličin, nelinearitami, dopravním zpožděním atd. Korektní zjištění vzájemných kauzálních vazeb (jejich orientovanost a nonanticipativnost, tj. důsledky zásadně nepředcházejí své příčiny [ZÍTEK, P. 1990]) mezi veličinami systému je pak úkolem velice náročným, rozsáhlým a dlouhodobým. V teorii a praxi pak dochází u jednotlivých profesí k následujícím zjednodušujícím přístupům: • technologové: nedostatečná znalost statistických a kybernetických metod • odborníci v oblasti ekonometrie: zaměření teorie i aplikací pouze na oblast ekonomie • matematici-statistici: teoretický přístup bez znalosti praktických problémů • technici zaměřeni na automatizaci: neznalost moderních metod vyvíjených v ekonometrii, znalost pouze základů matematické statistiky pro průřezová data, avšak dobrá znalost teorie náhodných procesů, metod a SW pro identifikaci soustav • všichni: malá spolupráce a poznání problémů druhých, nedostatečné spojení teorie a praxe. V příspěvku je ukázána užitečnost jednotlivých přístupů při zjišťování kauzality, statistické analýze a identifikaci lineárních dynamických soustav.
2 Popis lineárních dynamických soustav Obecně jsou dynamické soustavy technologických procesů (technologické agregáty a zařízení) nelineární a mnohorozměrové (MIMO - jsou charakterizovány mnoha fyzikálně-chemickými rovnicemi i veličinami a popsatelné pomocí soustavy obecně nelineárních diferenciálních/diferenčních a algebraických rovnic s rozloženými parametry a s dopravním zpožděním). Při analýze je snahou soustavy vhodně dekomponovat a chápat jako spojení dílčích lineárních (linearizovaných) a jednorozměrových subsoustav SISO/MISO. Budeme tedy uvažovat případ takovéto subsoustavy podle obr.1, přičemž uvažovaná struktura modelu je typu ARX (regresní model) [LJUNG, L. 1987], [TŮMA, J. 1998], [NOSKIEVIČ, P. 1999]: v
y0 u
S
y
+
y(M)
+ Obr. 1. Schéma struktury uvažovaného modelu (u – vstupní/akční veličina, S – analyzovaná regulovaná soustava, v – porucha, y – výstupní veličina soustavy, y0 – počáteční hodnota výstupní veličiny, y(M) – výstupní/regulovaná veličina modelu). Pro lineární dynamické soustavy se spojitými přenosy ve tvaru: Y ( s) (1) , G S ( s) = U (s) lze získat diskretizované modely (diskrétní přenosy invariantní vzhledem k přechodové funkci) pomocí vztahu [VÍTEČEK, A. 1988]: (2) G (s) G SC ( z ) = (1 − z −1 ) Z {L−1 { S } t = kT } , s kde t = k.T je diskrétní čas v násobcích k = {0, 1, 2, ...} periody vzorkování T. Diferenční rovnice diskretizovaných modelů integračních a proporcionálních soustav s dopravním zpožděním mají obecný tvar (obsahující nonanticipativnost veličin): n m (3) y t = ∑ α i yt − i + ∑ β j u t − j , i =1
j =0
kde α i , β i jsou koeficienty, i, j = k.T je diskrétní zpoždění v násobcích k = {0, 1, 2, ...} periody vzorkování T, n ≤ m jsou meze sumace (stupně polynomů diskrétního přenosu soustavy). Nezpožděný vstup ut se vyskytuje pouze u ideálních proporcionálních soustav, u reálných soustav je vždy zpožděný (předbíhá výstup). Pokud označíme odhady koeficientů: ~ ~ (4) ai = α i , b j = β j , bude mít regresní diferenční rovnice pro uvažovanou strukturu modelu tvar: n m (5) (M ) yt = ∑ ai yt −i + ∑ b j ut − j + vt + et , i =1
j =0
kde et je chyba regresního odhadu (v uvažovaném případě zahrnující i numerickou chybu způsobenou určitou velikostí kroku a metodou integrace při generování odezvy soustavy).
3 Model jednoduché lineární dynamické soustavy Sp1 Srovnání přístupů analýzy bude dokumentováno na příkladu jednoduché lineární dynamické soustavy S typu Sp1 (Soustava proporcionální se setrvačností 1. řádu) se spojitým přenosem (6) k1 G S (s) = , T1 s + 1 na níž působí pouze vstupní náhodný signál u a aditivní porucha výstupu v je nulová. Uvažovaná soustava Sp1 má diskrétní přenos: T (7) − k1 (1 − c) z −1 T1 G SC ( z ) = c e , = , 1 − cz −1 s odpovídající diferenční rovnicí (obsahující nonanticipativnost veličin): (8) yt = α1 yt −1 + β1ut −1 , přičemž pro koeficienty platí vztahy: T (9) − T1 α 1 = c = e , β 1 = k1 (1 − c) . Odhady parametrů soustavy z odhadů regresních koeficientů lze získat (viz i obr.2): ~ ~ (10) b T k1 = 1 , T1 = − . 1 − a1 ln(a1 ) 3 k1( a1 , 0.1)
15 T1( a1 , 0.5)
2
k1( a1 , 0.2)
T1( a1 , 1 )
k1( a1 , 0.3) 1
T1( a1 , 2 )
0
0.5 a1
1
10
5
0
0.5 a1
Obr. 2. Grafy funkcí přepočtů koeficientů: k1(a1, b1) a T1(a1,T) Simulace modelu byla uskutečněna v programu Matlab 4.2c s následujícími parametry: • parametry soustavy: k1 = 2, T1 = 10 s, y0 = 10
1
•
vstupní signál:
• • •
porucha: parametry simulace: metoda simulace:
u ~ N ( µ u , σ u ), µ u = y 0 / k1 = 5, σ u = 1 , seed(u) = 123456 v = ε ~ N (0, σ ε ), σ ε = 0, seed (ε ) = 321 perioda vzorkování T = 1 s, doba simulace tf = 5.T1 = 50 s lsim(A, B, C, D, u, t, x0), x0 = y0.
4 Srovnání různých přístupů analýzy a identifikace Srovnáme několik přístupů k řešení otázky existence závislosti mezi signály u a y(M) uvažovaného modelu s podpůrnými cíly: - zjištění orientovanosti a nonanticipativnosti kauzality (přičemž apriori orientovanost známe: příčina = u, následek = y(M), což ale v praxi nebývá zdaleka samozřejmé), - kvalifikace a kvantifikace závislosti mezi signály - parametrická identifikace soustavy, tj. zjištění hodnot parametrů soustavy (strukturální identifikaci není třeba provádět, jelikož struktura modelu soustavy je známá - což ale v praxi nebývá také zdaleka samozřejmé). Jsou to přístupy: • expertní – tj. „okometrické” (vizuální) posouzení závislosti na základě grafického průběhu veličin („podívám se a vidím!?”) • klasický statistický – použití regresní analýzy pro průřezová data • teorie náhodných procesů a statistické analýzy časových řad • ekonometrie a moderní statistické analýzy časových řad • identifikace systému odhadem parametrů modelu. 4.1 Expertní (vizuální) přístup Na základě grafického znázornění (obr. 3) uvažovaných stacionárních veličin (které jsou při ustáleném chodu kontinuálních technologických procesů nejčastější) se pozorovateli zdá, že mezi veličinami závislost (spíše) neexistuje: ) *+, ( ' % #
!
"!
#!
$!
%!
&!
$!
%!
&!
/ 0 *+, " ! .& "! - .& -
!
"!
#!
+123 4
Obr. 3. Časové průběhy vstupní a výstupní veličiny modelu (Matlab) Pro lepší posouzení (ne)závislosti a potvrzení hypotézy si expert ještě může zobrazit korelační diagram (bodový graf) veličin u a yM (obr. 4), který mu hypotézu nezávislosti „potvrdí” (díky skutečnosti že graf má tvar “mraku” bodů charakteristického pro nezávislé veličiny).
Bodový graf veličin yM - u 10.6 10.4 10.2
yM
10.0 9.8 9.6 9.4 9.2 2.5
3.5
4.5
5.5
6.5
7.5
8.5
u
Obr. 4. Korelační diagram výstupní (yM) a vstupní (u) veličiny modelu (STATISTICA) Závěr: Jde o jev (problém) tzv. zdánlivé nezávislosti (nekorelovanosti, nonkorelace) stacionárních veličin, kterému může expert (technolog, technik) lehce podlehnout. Duálním problémem je naopak zdánlivá závislost (korelace), která se vyskytuje hlavně u nestacionárních veličin a má obdobné neblahé důsledky při subjektivním vizuálním hodnocení. 4.2 Základní klasický statistický přístup Klasický statistický přístup vychází z použití regresní analýzy pro průřezová data - tak jak je doposud vyučován na univerzitách [ANDĚL, J. 1985] (s výjimkou výuky specialistů) a v praxi běžně používán. Přitom je možné použít regresní rovnici s absolutním členem a bez něj. a) Regresní rovnice s absolutním členem Regresní rovnice (model) bude mít v tomto případě zjednodušený (statický) tvar: (M ) yt = a 0 + b0 u t + et = a + bu t + et .
(11)
A. Výsledky lineární regrese z běžně používaného software (tabulkové procesory – např. Excel, základní statistické programy – např. Statgraphics) jsou uvedeny v tab. 1: Tab. 1. Základní výsledky jednoduché lineární regrese (Statgraphics) Parametr / statistika Hodnota Poznámka / závěr absolutní člen a 9.562 statisticky významný statisticky nevýznamný lineární člen b (0.053) statisticky nevýznamný Fisherův F-test modelu (2.356) 2 statisticky nevýznamný koeficient determinace R (0.046) Pozn.: Uvažovaná hladina významnosti α = 0,05, hodnoty v závorkách jsou na této hladině statisticky nevýznamné, celkový F-test modelu je v tomto případě ekvivalentní testu koeficientu determinace, výsledky byly zaokrouhleny na 3 desetinná místa. Závěr: Uvedený model lineární regrese a lineární člen jsou statisticky nevýznamné, tj. mezi veličinami neexistuje lineární závislost. Expert si tedy potvrdí svou domněnku zdánlivé
nezávislosti a téměř s jistotou bude konstatovat skutečnou nezávislost obou veličin – počítač a renomovaný statistický SW mu ji přece potvrdil! B. Trochu zkušenější expert pomocí kvalitního statistického SW navíc prozkoumá i výsledky analýzy reziduí (ověření jejich autokorelace, normality a homoskedasticity) – tab. 2. Tab. 2. Analýza reziduí jednoduché lineární regrese (EasyReg, TSP) Vlastnost Test Hodnota Poznámka / závěr pozitivní autokorelace autokorelace Durbin – Watson 0.400 normalita Jarque – Berra 1.871 akceptována homoskedasticita Breusch – Pagan 0.147 akceptována -”LM 0.178 akceptována Pozn.: U většiny běžných statistických programů lze normalitu a homoskedasticitu (konstantnost rozptylu) reziduí posoudit pouze graficky. Kritická hodnota Durbin-Watsonova testu autokorelace pro jednu vysvětlující proměnnou (abs. člen se nepočítá), n = 51 hodnot a předpoklad pozitivní autokorelace je např. podle [ANDĚL, J. 1993] dL( α = 0.05) = 1,50. Závěr: Z výsledků lze konstatovat, že rezidua jsou významně pozitivně korelována. Jak uvádí statistická literatura, důvodů může být více: nesprávně zvolený regresní model, neuvažování časových posunů tzv. zpožděných proměnných, nedostatečný/nevhodný počet/výběr regresorů apod. Zde je tedy první signalizace skutečnosti, že model je nesprávný. C. Mimořádně znalý a zkušený uživatel kvalitního statistického SW si prohlédne celkové hodnocení tzv. regresního tripletu (data + model + metoda) – tab. 3: Tab. 3. Testování regresního tripletu (ADSTAT, QC Expert) Objekt Model
Rezidua
Vlastnost významnost korektnost heteroskedasticita normalita autokorelace -”trend
Test
Hodnota
Fisher-Snedocorův Scottovo kritérium Cook – Weisbergův Jarque – Berraův Waldův Durbin-Watsonův znaménkový
2.356 -0.998 0.136 1.871 113.4 0.400 -3.566
Poznámka / závěr není významný je korektní zamítnuta potvrzena existuje pozitivní vykazují
Závěr: Rezidua jsou autokorelována, model je nesprávně sestaven a proto je statisticky nevýznamný. Autokorelace reziduí je způsobena proměnnou, která je taky autokorelována. Po ověření základních předpokladů statistického zpracování a tím i regresní analýzy (nezávislost, normalita, homogenita, vlivné body výběrů/regresorů) dostaneme výsledky - tab. 4 (ADSTAT, QC Expert): Tab. 4. Ověření základních předpokladů Veličina / u yM vlastnost normalita prokázána prokázána autokorelace nezávislost prokázána trend není není homogenita přijata přijata vlivné body nejsou nejsou
Z výsledků je zřejmé, že veličina yM je pozitivně autokorelovaná po řád 2 (obr. 5) a způsobuje problémy při klasické regresní analýze průřezových dat.
Obr. 5. Graf normované autokorelační funkce veličiny yM (QC Expert) Zde možnosti klasické regresní analýzy pro průřezová data prakticky končí, protože „autokorelovaná data se zahazují, nebo popíšou!” [MELOUN, M. 1995]: - většinou je snaha popsat chybu regresního odhadu pomocí modelu AR(1), což příliš neosvětlí řešení problému [WONNACOT, T. H. & WONNACOT, R. J. 1993], - dalšími možnostmi řešení je použít Aitkenovu zobecněnou MNČ a Cochran-Orcuttův iterativní postup [HEBÁK, P. & HUSTOPECKÝ, J. 1987], [GARAJ, V. & ŠUJAN, I. 1980], [GREŃ, J. 1984], ovšem postupy jsou pracné a získané výsledky problematické z hlediska jejich interpretace [MORÁVKA, J. & MRAJCA, V. 1999]. b) Regresní rovnice bez absolutního členu Další možností pro srovnání je použití lineárního modelu bez absolutního členu (což však běžného uživatele stěží napadne a navíc statistická literatura to z různých racionálních důvodů ani nedoporučuje – viz např. [VÍŠEK, J. Á. 1998]): (M ) (12) yt = b0 u t + et = bu t + et . Takovýto model má však řadu nectností: koeficient determinace R2 ztrácí svůj původní význam a v různých programech bývá pro tento případ jinak stanoven (tab. 5). Tab. 5. Celkové výsledky pro rovnici bez absolutního členu (ADSTAT, QC Expert, TSP) Objekt Parametr / vlastnost Koeficient / test Hodnota Hodnocení regresní koeficient b je významný 1.914 významnost Fisher-Snedocorův 1469 je významný Model koeficient determinace R2 0.967 / ± 57.9! (je významný) korektnost Scottovo kritérium 0.00 je korektní heteroskedasticita Cook – Weisbergův 2.096 zamítnuta homoskedasticita LM 0.502 akceptována normalita Jarque – Berraův 0.153 potvrzena Rezidua autokorelace Waldův 0.185 není -”Durbin-Watsonův 1.873 není trend znaménkový 0.070 není Závěr: Výsledky jsou podstatně jiné a lepší, než v předchozím případě. Model je statisticky významný, korektní a rezidua et mají charakter Gaussovského náhodného procesu (bílého šumu). Regresní koeficient b je dobrým odhadem koeficientu přenosu soustavy k1. Tento
statický model však neumožňuje podchytit dynamiku soustavy a tak nelze odhadnout její časovou konstantu. Závěr: Klasická regresní analýza průřezových dat dává (částečně) správné výsledky pouze pro nezávislé, stejně rozdělené i homogenní veličiny a při správné volbě modelu. Logickým důsledkem je tedy další zkoumání uvedených veličin jako časových řad – náhodných procesů. 4.3 Přístup teorie náhodných procesů a analýzy časových řad Teorie náhodných procesů a statistická analýza časových řad používá pro charakterizaci signálů statistické funkce vyššího řádu, než jakými jsou základní charakteristiky (průměr, medián, rozptyl, šikmost, špičatost atd.). Tyto „vyšší” funkce popisují průběh veličin v časové a frekvenční oblasti. Pro jeden signál (proces, časovou řadu) jsou používány: – (normovaná) autokorelační funkce (AKF, nebo ACF – Auto Correlation Function) – (normovaná) parciální autokorelační funkce (PAKF, nebo PACF – Partial Auto Correlation Function) [CIPRA, T. 1986], [ARLT, J. 1999] – (výkonová) spektrální hustota, autospektrum (VSH, nebo PSD – Power Spectral Density), periodogram (jako odhad spektrální hustoty) [BENEŠ, J. 1961], [CIPRA, T. 1986] Pro dva signály (pro testování jejich vzájemné závislosti a vazby) se používá: – (normovaná) vzájemná korelační funkce (VKF, nebo CCF – Cross Correlation Function) – vzájemná (výkonová) spektrální hustota, křížové spektrum (VVSH, nebo CPSD) – koherenční funkce [ANDĚL, J. 1975], [CIPRA, T. 1986], [TŮMA, J. 1997] – funkce zisku [CIPRA, T. 1986] – fázová funkce (fázové spektrum) [ANDĚL, J. 1975], [CIPRA, T. 1986]. A. Prozkoumáme nejprve charakter jednotlivých signálů pomocí grafického znázornění (diagramů) uvedených funkcí: "
korelogramy autokorelačních a parciálních autokorelačních funkcí (Statgraphics):
Závěr: Na základě klasifikace podle metodologie Box-Jenkinse vycházející z tvaru a identifikačních bodů diagramů ACF a PACF [CIPRA, T. 1986] lze konstatovat, že časová řada u(t) má charakter náhodného procesu typu bílý šum, zatímco řada yM(t) má charakter AR(2), tj. autoregresního procesu 2. řádu, který naznačuje výstup ze setrvačné soustavy. "
diagramy (výkonových) spektrálních hustot signálů u a yM (STATISTICA):
Spectral analysis: U No. of cases: 50
Spectral Density
Hamming weights:.0357 .2411 .4464 .2411 .0357 5
5
4
4
3
3
2
2
1
1
0 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.45
0 0.50
Frequency
Spectral analysis: YM No. of cases: 50
Spectral Density
Hamming weights:.0357 .2411 .4464 .2411 .0357 0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.0 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.45
0.0 0.50
Frequency
Závěr: Na základě průběhu autospekter [BENEŠ, J. 1961] lze konstatovat, že časová řada u(t) má charakter přibližně širokopásmového náhodného procesu typu bílý šum, zatímco řada yM(t) má charakter nízkofrekvenčního šumu, tj. naznačuje výstup z proporcionální nekmitavé (aperiodické) soustavy se setrvačností určitého řádu, či z dolní propusti (dolnopropustného filtru). Z korelačních funkcí a spekter lze v tomto případě vypozorovat orientaci kauzality – po průchodu signálu dynamickou soustavou je výstupní signál více (auto)korelovaný a má užší spektrum (což vyplývá z duality funkcí ACF/PSD). B. Jelikož nás hlavně zajímá vzájemný vztah signálů u(t) a yM(t), dostaneme grafické znázornění (diagramy) příslušných funkcí (např. z programu STATISTICA): "
diagram normované vzájemné korelační funkce (VKF / CCF):
Cross-Correlation Function First : YM Lagged: U Lag Corr. -7 -.105 -6 -.081 -5 -.173 -4 -.180 -3 -.160 -2 -.212 -1 -.222 0 .2142 1 .5482 2 .3946 3 .2980 4 .1459 5 .0352 6 .0168 7 .0993
S.E. .1508 .1491 .1474 .1459 .1443 .1429 .1414 .1400 .1414 .1429 .1443 .1459 .1474 .1491 .1508 -1.0
-0.5
0.0
0.5
1.0
Závěr: Z obrázku VKF je vidět skutečnost, že veličina u „předbíhá” veličinu yM o 1 krok (maximum VKF je pro hodnotu předstihu +1, ale statisticky významná je i hodnota VKF(+2)), neboli výstup yM je zpožděn za vstupem u o 1 krok. Opačná závislost není statisticky významná a tak orientace kauzality je jednoznačná: u → yM. V regresní rovnici pro yM(t) se proto musí objevit na pravé straně vstupní veličina předsunuta o 1 krok, tj. u(t-1). "
diagram vzájemné spektrální hustoty: Cospectral Density X:U Y:YM No. of cases: 50 (trunc.)
Cospectral Density
Bartlett weights:0.000 .2500 .5000 .2500 0.000 0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.0
0.0
-0.1 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.45
-0.1 0.50
Frequency
Závěr: Z diagramu křížového spektra vyplývá závislost signálů pouze pro nízké frekvence. "
koherenční diagram:
Squared Coherency X:U Y:YM No. of cases: 50 (trunc.)
Squared Coherency
Bartlett weights:0.000 .2500 .5000 .2500 0.000 1.0
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0.0 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.45
0.0 0.50
Frequency
Závěr: Z koherenčního diagramu vyplývá závislost signálů pro nižší a střední frekvence. Zde platí (koherenční funkce je blízká hodnotě 1), že sledovaný model je lineární, má jeden vstup a signály na vstupu/výstupu nejsou ovlivněny aditivním šumem [TŮMA, J. 1997]. diagramy funkce zisku (Gain): Gain of X over Y X:U Y:YM No. of cases: 50 (trunc.) Bartlett weights:0.000 .2500 .5000 .2500 0.000
Gain of X over Y
"
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.0 0.00
0.05
0.10
0.15
0.20
0.25 Frequency
0.30
0.35
0.40
0.45
0.0 0.50
Gain of Y over X X:U Y:YM No. of cases: 50 (trunc.)
Gain of Y over X
Bartlett weights:0.000 .2500 .5000 .2500 0.000 30
30
25
25
20
20
15
15
10
10
5
5
0 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.45
0 0.50
Frequency
Závěr: Z diagramů funkce zisku lze konstatovat, že hodnoty lineárních regresních koeficientů v obou variantách kmitočtového modelu bez absolutního členu jsou výrazně závislé na frekvenci [CIPRA, T. 1986]. Vliv veličiny u na veličinu yM je největší pro nízké frekvence a klesá s rostoucí frekvencí. Naopak vliv veličiny yM na veličinu u je největší pro nejvyšší frekvence a klesá s klesajícím kmitočtem. Pro nízké kmitočty je tento vliv zanedbatelný a z těchto skutečností můžeme usuzovat na orientaci kauzality: u → yM. "
fázový diagram (diagram fázového spektra): Phase Spectrum X:U Y:YM No. of cases: 50 (trunc.)
Phase Spectrum
Bartlett weights:0.000 .2500 .5000 .2500 0.000 2.0
2.0
1.5
1.5
1.0
1.0
0.5
0.5
0.0
0.0
-0.5 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.45
-0.5 0.50
Frequency
Závěr: Z obrázku fázového diagramu je zřejmé, že strmější průběh fázové funkce pro nízké frekvence odpovídá časovému zpoždění signálu yM za signálem u. Pro frekvence střední má fázová funkce mírný, téměř konstantní průběh a odpovídá úhlovému zpoždění signálu yM za signálem u. Uvedený průběh fázové funkce je typický pro ekonomické časové řady [CIPRA, T. 1986] a pro dynamické systémy (filtry) typu „dolní propust”. Pro vysoké frekvence jsou zpoždění vzhledem k průběhu koherenčního diagramu nevýznamné.
Závěr: Z uvedených diagramů lze odpozorovat pouze následující kvalitativní skutečnosti: • orientaci (směr) kauzality signálů: u → yM (VKF, funkce zisku, ACF/PACF a PSD) • nonanticipativnost: veličina u předbíhá veličinu yM o jeden krok (VKF) • charakter soustavy je typu „dolní propust” (PSD, fázová funkce) • soustava je lineární (koherenční funkce). Pro určení parametrů soustavy (kvantitativní analýza) je třeba použít jiných přístupů. 4.4 Ekonometrická a moderní statistická analýza časových řad Tento přístup předpokládá stanovení nebo znalost struktury tzv. dynamického regresního modelu se zpožděnými proměnnými. Řešením dříve odvozené dynamické regresní rovnice (8) prostřednictvím vhodného statistického SW získáme následující výsledky (tab. 6): Tab. 6. Výsledky regresního tripletu pro regresní model soustavy Sp1 (ADSTAT, TSP) Objekt Parametr / vlastnost Koeficient / test Hodnota Hodnocení a1 0.950 statisticky významný regresní koeficienty 0.101 statisticky významný b1 Model významnost Fisher-Snedocorův 250 je významný korektnost Scottovo kritérium -0.940 je korektní heteroskedasticita Cook – Weisbergův 0.374 zamítnuta homoskedasticita LM 0.488 akceptována normalita Jarque – Berraův 0.264 potvrzena Rezidua autokorelace Waldův 0.0004 není -”Durbin-Watsonův 1.989 není -”Durbinův-h 0.019 není trend znaménkový 0.765 není Pozn.: Pro testování autokorelace v regresním modelu se zpožděnými vysvětlovanými proměnnými není vhodné použít Durbin-Watsonův test, ale Durbinův-h test [ARLT, J. 1999]. Grafický průběh výstupu soustavy i modelu včetně průběhu reziduí et regresního odhadu je na následujícím obrázku (obr. 6): yM
model
10
9.5 5 3
10
15
20
25
30
35
40
45
50
10
15
20
25
30
35
40
45
50
rezidua
2 1 0 -1 -2 5
Obr. 6. Průběh výstupu soustavy i modelu včetně průběhu reziduí (PcGive)
Závěr: Model je po všech stránkách statisticky významný a korektní. Odhady parametrů soustavy (přepočtené z regresních koeficientů) mají hodnoty: ~ ~ (13) k1 ≈ 2.02, T1 ≈ 19.5 s , a byly stanoveny s relativními chybami asi +1% a +95 %. Po provedení analýzy orientace kauzality časových řad pomocí analýzy impuls-reakce (I-R) [ARLT, J. 1999] pro vektorový autoregresní model 1.řádu VAR(1) systému z(t)=[u(t), yM(t)]T, dostaneme grafické průběhy reakcí časových řad na impuls v řadě u(t) - obr.7, 8:
Obr. 7. Reakce řady u(t) na impuls v řadě u(t) (EasyReg)
Obr. 8. Reakce řady yM(t) na impuls v řadě u(t) (EasyReg) Závěr: Reakce řady yM(t) na impuls v yM(t) a reakce řady u(t) na impuls v řadě yM(t) jsou identicky nulové. Na obrázcích je vidět, že řada yM(t) reaguje se zpožděním jednoho kroku na impuls v řadě u(t), co odpovídá odezvě soustavy Sp1. V případě analýzy systému s opačným pořadím časových řad, tj. z(t) = [yM(t), u(t)]T však dostaneme nejasné výsledky. Pro správnou analýzu impuls-reakce je nutné dopředu stanovit pořadí procesů zařazených do systému VAR na základě věcné analýzy – což je nevýhodou při určování orientace kauzality ve srovnání s funkcí VKF a funkcí zisku, z jejichž diagramů lze orientaci kauzality určit bez předchozích znalostí. Závěr: Při znalosti struktury modelu soustavy jsou výsledky dynamických regresních rovnic naprosto významné a korektní. Při správně odhadnutém (odzkoušeném) směru kauzality lze odhadnout také zpoždění výstupu (důsledku) za vstupem (příčinou, akcí).
4.5 Identifikace systému odhadem parametrů modelu Pomocí přístupů a metod uvedených v [LJUNG, L. 1987], [TŮMA, J. 1998], [NOSKIEVIČ, P. 1999] lze v programovém systému Matlab/Simulink – System Identification Toolbox (nebo v dalších obdobných programech) získat pro strukturu ARX prostřednictvím příkazů: th_arx = arx([yM u], [na nb nk]); % kde na=1, nb=1, nk=1 present(th_arx) následující odhady koeficientů a kvality (přesnosti) získaného modelu: a1 = 0.9499, se(a1) = 0.0074 ⇒ a1/se(a1) = 128.4 b1 = 0.1008, se(b1) = 0.0144 ⇒ b1/se(b1) = 7.0 Loss function = 0.008772, FPE = 0.009488. Pozn.: Koeficient a1 je uveden s opačným znaménkem s ohledem na srovnání se statistickým výpočtem (struktury obou rovnic jsou jinak sestaveny). Označení se(p) znamená směrodatnou chybu (standard error) odhadu koeficientu p. Poměr se(p)/p má význam pro hodnocení statistické významnosti koeficientů – přibližně platí, že poměr > 2 znamená statistickou významnost koeficientu (pro přesné ohodnocení by bylo třeba použít kritických hodnot Studentova t-rozdělení). Odhadnuté koeficienty jsou tedy statisticky vysoce významné. Loss function znamená ztrátovou funkci a FPE (final prediction error) Akaikeho konečnou chybu predikce. Oboje mají význam pro porovnání a hodnocení přesnosti modelu – hodnoty obou kritérií byly minimální pro metodu arx v porovnání s dalšími metodami armax, oe, iv4 a bj (seřazeno vzestupně). Závěr: Při dopředu neznámých typech soustav a poruch lze prostřednictvím uvedeného SW určit pro model arx pomocí příkazů arxstruc, selstruc, struc nejvhodnější strukturu charakterizovanou počtem kroků dopravního zpoždění a stupni polynomů koeficientů [TŮMA, J. 1998]. Pro hodnocení statistické významnosti získaných koeficientů je třeba stanovit poměry parametrů a jejich směrodatných chyb. Parametry spojitých soustav lze získat přepočtem a srovnáním příslušného spojitého a diskrétního přenosu. Pro náš případ jsou jejich hodnoty: ~ ~ (14) k1 ≈ 2.02, T1 ≈ 19.5 s , přičemž parametry byly stanoveny s relativními chybami asi +1% a +95 %. Odhady parametrů soustavy jsou identické s odhady získanými pomocí ekonometrické analýzy.
5 Závěr Analýza a identifikace dynamických soustav v průmyslové praxi není jednoduchou záležitostí. Hodnocení informací o soustavách na základě měřených veličin naráží na mnohá omezení a problémy. Nesprávné výsledky (i jejich interpretace) a selhání přístupů se objevuje u: • expertního odhadu v případě zdánlivé nonkorelace (a zdánlivé korelace) veličin • klasických statistických metod v případě nesprávně zvoleného regresního modelu a/nebo nedostatečných znalostí experta (pracujícího navíc s nevhodným SW). Správných výsledků s ohledem na definované cíle, tj. zjištění kauzality (aspektů orientovanosti a nonanticipativnosti), kvalitativní (charakter soustavy, linearita) i kvantitativní analýza (odhady regresních koeficientů) a identifikace (strukturální a parametrická) dynamických soustav lze dosáhnout pouze vhodnou kombinací výše popisovaných přístupů – tab. 7.
Tab. 7. Shrnutí výsledků popisovaných přístupů Kauzalita Výsledky, Přístup Poznámka závěry orientace nonanticipace Expertní nesprávné model s abs. autokorelace nesprávné členem reziduí Klasický statistický částečně neúplné – bez abs. členu správné statický model Teorie náhodných procesů a + (CCF, pouze +dolní propust, + (CCF) analýza časových řad Gain) kvalitativní linearita Ekonometrie a moderní nutná znalost (+) (I-R) (+) (I-R) správné metody analýzy časových řad struktury pro ARX Identifikace soustav odhadem + správné i zjištění parametrů modelu struktury
6 Literatura ANDĚL, J. 1975. Statistická analýza časových řad. 1.vyd. Praha : TKI SNTL, 1975. 282 s. ANDĚL, J. 1985. Matematická statistika. 2.vyd. Praha : SNTL/ALFA, 1985. 352 s. ANDĚL, J. 1993. Statistické metody. 1.vyd. Praha : Matfyzpress MFF UK Praha, 1993. 246 s. ARLT, J. 1999. Moderní metody modelování ekonomických časových řad. 1.vyd. Praha: Grada Publishing, s.r.o., 1999. 312 s. ISBN 80-7169-539-4. BENEŠ, J. 1961. Statistická dynamika regulačních obvodů. 1.vyd. Praha : SNTL, 1961. 336 s. BIERENS, H. 1999a. Cointegration Analysis [online]. Pensylvania State University, PA, 1999.
29 s. BIERENS, H. 1999b. Nonparametric Nonlinear Co-Trending Analysis, With an Application to Interest and Inflation in the U.S. [online]. Pensylvania State University, PA, 1999.
39 s. BIERENS, H. 1999c. Free Econometrics Software for Easy Regression Analysis [online]. Pensylvania State University, PA, 1999.
CIPRA, T. 1986. Analýza časových řad s aplikacemi v ekonomii. 1.vyd. Praha : SNTL/ALFA, 1986. 248 s. GARAJ, V. & ŠUJAN, I. 1980. Ekonometria. 1.vyd. Bratislava : ALFA/SNTL, 1980. 288 s. GREŃ, J. 1984. Statystyka matematyczna – modele i zadania. VIII.vyd. Warszawa : PWN, 1984. 326 s. ISBN 83-01-03699-0. HEBÁK, P. & HUSTOPECKÝ, J. 1987. Vícerozměrné statistické metody s aplikacemi. 1.vyd. Praha : SNTL/ALFA, 1987. 456 s. ISBN 80-01-01076-7. HUŠEK, R. & WALTER, J. 1976. Ekonometrie. 1.vyd. Praha : SNTL, 1976. 264 s. JOHANSEN, S. 1997. Mathematical and Statistical Modelling of Cointegration. In Sborník EUI (European University Institute) Florence, Italy : Working Paper ECO No. 97/14, 1997, 16 s. LJUNG, L. 1987. System Identification : Theory for the User. 1.vyd. Englewood Cliffs New Jersey : PTR Prentice Hall, 1987. 519 s. ISBN 0-13-881640-9. MELOUN, M. & MILITKÝ, J. 1994. Statistické zpracování experimentálních dat. 1.vyd. Praha : PLUS, 1994. 839 s. ISBN 80-85297-56-6. MELOUN, M. 1995. Postupy regresní diagnostiky. Přednáška na semináři Analýza dat ‘95/II, Lázně Bohdaneč, 21.-24.11.1995. MORÁVKA, J. & MRAJCA, V. 1999. Ekologická optimalizace provozu A2 SP4. (Analytická studie projektu č. 5098030). Třinec : TŽi a.s. SA - Středisko automatizace, červen-říjen 1999. 40 s.
NOSKIEVIČ, P. 1999. Modelování a identifikace systémů. 1.vyd. Ostrava : MONTANEX, 1999. 276 s. ISBN 80-7225-030-2. TŮMA, J. 1997. Zpracování signálů získaných z mechanických systémů užitím FFT. 1.vyd. Praha : nakladatelství Sdělovací technika, 1997. 174 s. TŮMA, J. 1998. Složité systémy řízení I. – Regulace soustav s náhodnými poruchami. 1.vyd. Ostrava : skripta FS VŠB-TU Ostrava, 1998. 158 s. VÍTEČEK, A. 1988. Matematické metody automatického řízení (Transformace L a Z). 1.vyd. Ostrava : skripta FSE VŠB Ostrava, 1988. 156 s. VÍTEČEK, A. & WAWRZICZKOVÁ, M. 1988. Teorie automatického řízení. 1.vyd. Ostrava : skripta PGS VŠB Ostrava, 1988. 175 s. VÍŠEK, J. Á. 1998. Statistická analýza dat. I.vyd. Praha : skripta FJFI ČVUT- Praha, 1998. 187 s. ISBN 80-01-01735-4. WONNACOT, T. H. & WONNACOT, R. J. 1993. Statistika pro obchod a hospodářství. 1.vyd. Praha : Victoria Publishing, 1993, 891 s. ISBN 80-85 605-09-0. ZÍTEK, P. 1990. Simulace dynamických systémů. 1. vyd. Praha : SNTL, 1990. 420 s. ISBN 8003-00330-X.
