UNIVERZITA PALACKÉHO V OLOMOUCI PŘÍRODOVĚDECKÁ FAKULTA KATEDRA MATEMATICKÉ ANALÝZY A APLIKACÍ MATEMATIKY
BAKALÁŘSKÁ PRÁCE Design experimentu
Vedoucí diplomové práce: Mgr. Jaroslav Marek, Ph.D. Rok odevzdání: 2010
Vypracovala: Kristina Axmanová ME, III. ročník
Prohlášení Prohlašuji, že jsem diplomovou práci zpracovala samostatně pod vedením pana Mgr. Jaroslava Marka, Ph.D. s použitím uvedené literatury.
V Olomouci dne 14. dubna 2010
Poděkování Ráda bych poděkovala svému vedoucímu bakalářské práce Mgr. Jaroslavu Markovi, Ph.D., za cenné rady, odborné informace a čas strávený při konzultacích. Dík také patří rodině a přátelům, kteří mi pomáhali s experimentem a podporovali mě po celou dobu studia.
Obsah Úvod
4
1 Plánování experimentů 1.1 Historie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Základní pojmy . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Příprava experimentu . . . . . . . . . . . . . . . . . . . . . . . . .
5 5 6 7
2 Typy faktorových plánů 2.1 Úplný faktorový experiment . . . . 2.2 Jednofaktorový experiment . . . . . 2.3 Částečný faktorový experiment . . 2.3.1 Poloviční plán . . . . . . . . 2.3.2 Saturované (nasycené) plány 2.3.3 Středové plány . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
9 9 10 11 12 14 15
3 Průběh experimentu 3.1 Úplný faktorový experiment s výsledky měřění 3.1.1 Výpočet efektu faktorů . . . . . . . . . 3.1.2 Test významnosti efektu . . . . . . . . 3.2 Poloviční plán . . . . . . . . . . . . . . . . . . 3.2.1 První polovina plánu . . . . . . . . . . 3.2.2 Druhá polovina plánu . . . . . . . . . . 3.2.3 Shrnutí . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
16 16 17 19 22 22 24 25
4 Regresní model 4.1 Regresní model Box–Behnken design 4.1.1 Vzorový příklad . . . . . . . . 4.1.2 Vlastní příklad . . . . . . . . 4.2 Dodatečná měření . . . . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
25 26 26 29 31
. . . .
. . . .
. . . .
. . . .
. . . .
5 Závěr
35
Literatura
36
Úvod V praxi se často setkáváme s úlohami, které vedou k regresní analýze. Problémem ale může být najít vhodnou funkci popisující vztah mezi jednou vysvětlovanou proměnnou a vysvětlujícími proměnnými. K nalezení vhodného regresního modelu nám slouží několik metod, mezi které patří i metoda plánování experimentu označovaná jako DOE (Design of Experiment). Účelem plánování experimentu je zkoumat vliv faktorů na výstupní proměnnou a stanovit, které z nich, resp. které jejich interakce, jsou statisticky významné a tedy ovlivňují sledovanou závislou proměnnou. Pro velký počet vysvětlujících proměnných byly zkonstruovány různé faktorové plány, pomocí nichž můžeme zredukovat počet zkoumaných faktorů. Použití těchto faktorových plánů umožňuje výběr vhodného regresního modelu a nalezení kvalitnějších odhadů neznámých parametrů aproximující funkce a k celkově efektivnějšímu způsobu získání výsledků experimentu. Jejich aplikace pak může vést ke snížení nákladů na výrobu. a k celkově efektivnějšímu způsobu získání výsledků experimentu. Prvním cílem práce je prostudovat vybrané partie z teorie plánování experimentu a zejména různé typy faktorových plánů. Druhým cílem je demonstrování jejich uplatnění na reálné úloze, která bude věnována „kvalitě výroby papírových vrtulníkůÿ. Sledovaným znakem kvality bude doba letu. Užitím faktorových plánů určíme významnosti jednotlivých faktorů (šířka, délka papíru a úhel startu) a jejich interakcí. Po nalezení regesního modelu najdeme parametry výroby papírového vrtulníku, které nám prodlouží dobu jeho letu.
4
1. Plánování experimentů 1.1. Historie Teorie plánování experimentu je oblast matematické statistiky, která vznikla na počátku 20. století. Text k této partii vznikl kompilací textů z literatury [3, 6, 8]. Plánování experimentu je spojeno s vynikajícím anglickým matematikem Ronaldem Fisherem. Od roku 1919 zpracovával výsledky experimentální práce v Rothamsteadské královské zemědělské výzkumné stanici, které vznikaly po dobu šedesáti let. Při tomto výzkumu položil teoretické základy analýzy rozptylu a navrhování experimentu. Sedmiletá práce byla završena sepsáním jeho první knihy s názvem Statistické metody pro výzkumné pracovníky, která byla přeložena do několika jazyků. Další velmi významnou knihou, napsanou v roce 1935, je Design of Experiments. S novým experimentálním návrhem, při kterém je počet pokusů násobkem čtyř, přicházejí v roce 1940 vědci Plackett a Burman. V 70. letech se mezi vědci objevuje statistik George Box, který navrhl několik modelů experimentu, např. Box-Jenkins models, Box-Cox transformace, BoxBehnken designs. Douglas C. Montgomery, působící na státní univerzitě v Arizoně, napsal v roce 1976 učebnici Design and Analysis of Experiments, obsahující základy navrhování experimentu. Tuto publikaci dodnes využívají nejen praktici, ale i tvůrci různých softwarů. V posledních letech se požaduje zvyšování jakosti výrobků a procesů, kterých se docílí propracovaným systémem metod, např. Taguchiho metoda. Japonský odborník Genichi Taguchi přistupuje k navrhování experimentu pro veřejnost srozumitelnou formou, kdy není vyžadována hluboká znalost statistiky. V praxi se tato metoda používá ve velké míře zejména v zahraničí, nejvíce pak v Japonsku a o něco později v USA. U nás se s ní setkáváme v devadesátých letech v automobilovém průmyslu. 5
Z popsaného vývoje v oblasti designu experimentu je zřejmé, že bylo navrženo obrovské množství různých metod. Některé z nich si v této práci popíšeme a vyzkoušíme v praxi.
1.2. Základní pojmy Experiment je soustava pokusů, ve kterých měníme obvyklé pracovní podmínky s cílem nalézt nejlepší pracovní postupy a také získat poznatky o vlastnostech výrobku a výrobního procesu, viz [7]. Experimentální postupy rozdělujeme na dva druhy: a) experimenty neplánované, b) experimenty plánované. Plánované experimenty se řídí určitým plánem experimentu. Plán experimentu stanovuje počet pokusů, ze kterých se experiment skládá, podmínky, za kterých se jednotlivé pokusy uskuteční a pořadí pokusů. Rozlišujeme význam pojmů pokus a experiment. Pokus je zjištění hodnoty ukazatele kvality za určitých, předem plánovaných podmínek výroby a experiment je systém všech pokusů. Upřesníme pojem nejlepší pracovní postup. Označíme-li sledovaný ukazatel kvality Y a faktory, které jej ovlivňují A, B, C, D, . . . se mohou pohybovat na různých úrovních (pro faktor A na úrovních A1 , A2 , A3 , . . . , pro faktor B na úrovních B1 , B2 , B3 , atd.), pak cílem plánování experimentu je a) rozhodnout, které z faktorů A, B, C, D, . . . významným způsobem ovlivňují ukazatel kvality Y , b) určit optimální úrovně významných faktorů (viz [7]). Cílem plánování experimentu je určit, zda určité vstupní proměnné (faktory) mají vliv na sledovanou proměnnou, viz [4]. • Faktory — jsou vstupní proměnné, které ovlivňují hodnoty sledované veličiny. Pro každý faktor jsou stanoveny úrovně (minimálně dvě), které vedou k dosažení optima sledované veličiny. 6
• Sledovaná veličina — pomocí této veličiny vyjádříme výsledek experimentu. Různé výsledky získáme změnou úrovní faktorů a jejich interakcí. Základní principy, které musí být zachovány při plánování experimentu: • replikace - je opakování experimentu při stejné úrovni faktorů. Opakováním odhalíme případné nepřesnosti měření a tím můžeme zvýšit spolehlivost závěru. • náhodnost - náhodné proměnné mohou ovlivňovat výsledek experimentu, a proto je nutné jednotlivé pokusy provést v náhodném pořadí. • rozdělení do bloků - experimentální jednotky rozdělíme do bloků, které vykazují podobné vlastnosti. Tím částečně zamezíme vlivu zdrojů variability. Podle typu zkoumaného problému můžeme stanovit jeden z následujících typů experimentu: • jednofaktorový experiment, • úplný faktorový experiment, • částečný faktorový experiment.
1.3. Příprava experimentu V knize [1] na str. 71 je uveden příklad, ve kterém se hledají mezi délkou pružiny, tloušťkou drátu a typem materiálu faktory resp. jejich interakce ovlivňující životnost pružiny. V této práci se budu věnovat obdobné úloze. Příklad 1. Budeme sledovat dobu letu vrtulníku na zem z výšky 20 m vyrobeného z archu papíru v závislosti na těchto faktorech: šířka papíru S, délka papíru D a startovací úhel vrtulníku U , viz obr. 1.
7
Obrázek č. 1: Výroba vrtulníku (viz [11]) Řešení: Určíme si maximální a minimální úrovně faktorů, které zapíšeme do tabulky.
faktor
Označení
Šířka papíru Délka papíru Úhel spuštění vrtulníku
Dolní úroveň Horní úroveň − 8 13 −90◦
S D U
+ 11 16 90◦
Tabulka č. 1: Faktory a jejich uvažované úrovně Plán experimentu lze sestavit více způsoby. Některé z nich si ukážeme v následující části. 8
2. Typy faktorových plánů 2.1. Úplný faktorový experiment Mezi nejpoužívanější plány experimentu patří úplný faktorový plán, který v daném případě vypadá takto:
pokus
S
D
U
1 2 3 4 5 6 7 8
8 11 8 11 8 11 8 11
13 13 16 16 13 13 16 16
90◦ 90◦ 90◦ 90◦ −90◦ −90◦ −90◦ −90◦
Y
Tabulka č. 2: Plán experimentu Y je výsledek pokusu. Plán experimentu je výhodnější psát pomocí takové symboliky, že dolní úroveň bude značena -1 a horní úroveň +1.
pokus
S
D
U
1 2 3 4 5 6 7 8
−1 1 −1 1 −1 1 −1 1
−1 −1 1 1 −1 −1 1 1
−1 −1 −1 −1 1 1 1 1
Y
Tabulka č. 3: Plán experimentu v kódovaných proměnných
9
Přepočet původních proměnných na tzv. kódované proměnné se provede takto:
xc =
xmax +xmin 2 xmax −xmin 2
x0 −
,
kde x0 = proměnná v původních jednotkách, xc = kódovaná proměnná, xmax = horní úroveň x, xmin = dolní úroveň x. Tento přepočet můžeme provést nejen pro krajní hodnoty xmax (= +1) a xmin (= −1). Počet pokusů úplného experimentu se vypočítá při k faktorech pomocí vztahu n = 2k . Při k = 3 faktorech je počet pokusů n = 23 = 8. Poznámka 1. Mohlo by se zdát, že úplný faktorový design je možné použít vždy. Ovšem je třeba si uvědomit, jaký je počet měření v závislosti na počtu faktorů. Pokud uvažujeme pouze dolní a horní úroveň každého z k faktorů, pak experimentů v úplném faktorovém designu je 2k . Například pro k = 5 potřebujeme 32 měření, pro k = 8 je dokonce nutných 256 měření. Je třeba si uvědomit, že ne vždy lze z důvodů ekonomických, časových či jiných tento s počtem faktorů narůstající počet měření uskutečnit. Tato situace vedla k nutnosti navrhnout postupy jiné a efektivnější. Některé z nich uvedeme v následující části. Situace je ještě horší, pokud spolu s dolní a horní úrovní sledovaných faktorů budeme ještě uvažovat střední úroveň, která je jejich průměrem a značí se symbolem „0ÿ. Potom je počet experimentů v úplném faktorovém modelu 3k , např. pro k = 5 je to 243 experimentů a pro k = 8 jde o 6561 měření.
2.2. Jednofaktorový experiment Už z názvu vyplývá, že se jedná o faktorový plán, ve kterém budeme zkoumat změnu jednoho faktoru. U zkoumaného faktoru se mění úroveň z dolní na horní, 10
přičemž u ostatních se drží na střední úrovni. Střední úroveň je vlastně průměrem dolní a horní úrovně, značíme 0. V našem příkladě, kdy zkoumáme tři faktory, by plán vypadal takto:
pokus
S
D
U
1 2 3 4 5 6
1 −1 0 0 0 0
0 0 1 −1 0 0
0 0 0 0 1 −1
Tabulka č. 4: Jednofaktorový plán
Počet pokusů jednofaktorového plánu se vypočítá při k faktorech pomocí vztahu n = 2k. Při k = 3 faktorech je počet pokusů n = 2 · 3 = 6.
2.3. Částečný faktorový experiment Na rozdíl od úplného faktorového experimentu, kde se sestavuje plán pro každý faktor, se u částečného faktorového plánu sestavuje plán jen pro několik faktorů (tzv. hlavní faktory) a ostatní faktory (tzv. vedlejší) se vyjádří jako kombinace hlavních faktorů. Částečné faktorové plány rozdělíme na: a) poloviční plány – plány s nejnižším stupněm snížení, b) Saturované (nasycené) plány – plány s nejvyšším stupněm snížení, c) středové plány – plány mezi nejnižším a nejvyšším stupněm snížení. Mějme úplný experiment 2k , kde 2 je počet úrovní faktorů a k je počet faktorů, pak 2k−p je částečný faktorový experiment, kde p je tzv. stupeň snížení, viz [7]. Nejmenší možné snížení počtu pokusů je snížení na polovinu, kde stupeň snížení p = 1, tj. 11
2k = 2k−1 , 2 kde 2k−1 je počet pokusů částečného faktorového experimentu. Např.: V plánu 26 máme n = 64 pokusů. Snížením o p = 1 se dostaneme na polovinu 26 = 26−1 = 25 = 32. 2 Stupeň snížení může být i p > 1. V takovém případě ale musí platit podmínka, že počet pokusů nesmí být menší než počet faktorů, tedy n ≥ k. V uvedeném případě, kdy k = 6 je nejvyšší možné snížení p = 3, tj. 26−3 = 23 = 8. Mezi nejnižším a nejvyšším snížením počtu pokusů může být ještě řada možností snížení, konkrétně pro k = 6 můžeme ještě snižovat na 26−2 = 24 = 16. 2.3.1. Poloviční plán Mějme 5 faktorů A, B, C, D, E. Dále označme jednotkový faktor I obsahující jen +. Smyslem polovičního plánu je určit hlavní faktory (např. A, B, C, D a vedlejší faktor E, který se určí jako kombinace hlavních faktorů, tj.
E = ABCDE. Faktor E se nazývá generátor plánu. Kombinace faktorů je slovo, které se skládá z písmen (faktorů).
12
Zavedeme vlastnosti operací s faktory: A · A = I, A · I = I · A = A, (A · B)C = A(B · C), A · B = B · A. Definiční rovnici určíme pomocí generátoru plánu s využitím uvedených vlastností: · E,
E = ABCD
E · E = E · ABCD, I = ABCDE. Pomocí této rovnice lze nalézt zaměnitelné interakce, což jsou interakce faktorů, které mají stejné posloupnosti znamének. Pokud např. chceme zjistit zaměnitelnou interakci k AB, vynásobíme definiční rovnici touto interakcí · AB,
I = ABCDE
AB · I = AB · ABCDE, AB = AA · BB · CDE, AB = CDE. Faktor E je možné sestavit mnoha způsoby, jelikož není nutné jej vyjadřovat pomocí všech faktorů. Např. a) E = CD, b) E = BCD a k nim příslušné definiční rovnice a) I = CDE, b) I = BCDE. Každá definiční rovnice má tzv. řešení plánu, které se zapisuje k typu plánu římskou číslicí jako index. Je to nejkratší slovo v definičních rovnicích. Řešením plánu pro a) a b) je tedy a) 25−1 III ,
b) 25−1 IV ,
kde římská číslice určuje počet písmen (faktorů) ve slově. 13
2.3.2. Saturované (nasycené) plány Jsou to plány s nejvyšším stupněm snížení, pro které platí, že počet vedlejších faktorů je stejný jako počet všech možných interakcí hlavních faktorů, což je splněno právě tehdy, když je počet faktorů 2r − 1, kde r = 2, 3, 4, .... Z toho vyplývá, že saturované plány nemůžeme použít pro každý počet faktorů. Např. Pro 6 faktorů zvolíme hlavní faktory (A, B, C) a vedlejší faktory (D, E, F ). Interakce hlavních faktorů jsou AB, AC, BC, ABC. Vidíme, že počet faktorů neodpovídá počtu interakcí, takže pro 6 faktorů nelze sestavit saturovaný plán. Uvažujme 7 faktorů A, B, C, D, E, F , G. Nejvyšší možné snížení je o p = 4, tj. 27−4 . Zvolíme hlavní faktory A, B, C a vedlejší faktory D, E, F , G, ze kterých pomocí interakcí AB, BC, AC, ABC vytvoříme generátory plánu (jejich počet je p = 4) například takto:
D = AB, E = BC, F = AC, G = ABC. Určíme definiční rovnice pro • vedlejší faktory D, E, F, G,
D = AB
·D
⇒
I = ABD,
E = BC
·E
⇒
I = BCE,
F = AC
·F
⇒
I = ACF,
G = ABC
·G
⇒
I = ABCG.
• součin dvou vedlejších faktorů DE, DF, DG, EF, EG, FG, D · E = AB · BC
⇒
DE = AC
· DE
⇒
I = ACDE,
D · F = AB · AC
⇒
DF = BC
· DF
⇒
I = BCDF,
D · G = AB · ABC E · F = BC · AC
⇒ ⇒
DG = C EF = AB 14
· DG · EF
⇒ ⇒
I = CDG, I = ABEF,
E · G = BC · ABC
⇒
EG = A
· EG
⇒
I = AEG,
F · G = AC · ABC
⇒
FG = B
· FG
⇒
I = BF G.
• součin tří vedlejších faktorů DEF, DEG, DFG, EFG, D · E · F = AB · BC · AC
⇒
DEF = I
D·E·G = AB·BC·ABC
⇒
DEG = B
·DEG
⇒
I = BDEG,
D·F ·G = AB·AC ·ABC
⇒
DF G = A
·DF G
⇒
I = ADF G,
E·F ·G = BC ·AC ·ABC
⇒
EF G = C
·EF G
⇒
I = CEF G.
DEF G = ABC
· DEF G
· DEF
⇒
I = DEF,
• součin všech vedlejších faktorů DEFG, D · E · F · G = AB · BC · AC · ABC
⇒
⇒ I = ABCDEF G. Počet definičních rovnic je 2p − 1 = 24 − 1 = 15. Pomocí nejvyššího stupně snížení jsme původních 27 = 128 pokusů snížili na 27−4 = 8 pokusů, takže počet zaměnitelných interakcí a faktorů je 128/8 = 16. Např. pro A je to A + BD + ABCE + CF + BCG + CDE + ABCDF + ACDG + BEF + + EG + ABF G + ADEF + ABDEG + DF G + ACEF G + BCDEF G, kde + znamená, že tyto kombinace mají stejnou posloupnost znamének. Z těchto dvojic má smysl ponechat jen interakce dvou: A + BD + CF + EG. 2.3.3. Středové plány Středové plány by měly být dobrou kombinací polovičních a saturovaných plánů. Řekli jsme si, že středových plánů může být více, a proto mezi nimi hledáme ten, který má méně pokusů než poloviční a obsahuje nejvýhodnější skupiny 15
zaměnitelných faktorů. Nejlepší zaměnitelnou interakcí je interakce obsahující největší počet faktorů. Skupina zaměnitelných faktorů závisí na volbě generátorů plánu. Pro faktory A, B, C, D, E, F , G porovnáme dvě různé volby v plánu 27−3 : a) generátory E=ABD, F=ABCD, G=ACD, pak definiční rovnice jsou I = ABDE = ABCDF = ACDG = CEF = BCEG = BF G = ADEF G 7−3 a stupněm řešení plánu je nejkratší slovo, tj. 2III .
b) generátory E=ACD, F=BCD, G=ABD, pak definiční rovnice jsou I = ACDE = BCDF = ABDG = ABEF = BCEG = ACF G = DEF G, stupněm řešení je 27−3 IV .
3. Průběh experimentu 3.1. Úplný faktorový experiment s výsledky měřění Z archu papíru o rozměrech S = 8; 11 a D = 13; 16 jsme vyrobili vrtulníky, které jsme pouštěli z 20 m na zem pod úhlem U = −90◦ ; 90◦ . Každý pokus jsme opakovali dvakrát a naměřené výsledky jsou zaznamenány v Tab. 4, kde Y vyjadřuje průměr obou měření.
pokus
S
D
U
Y1
Y2
Y
1 2 3 4 5 6 7 8
− + − + − + − +
− − + + − − + +
− − − − + + + +
15,0 11,0 15,5 14,6 14,4 14,5 15,8 13,8
14,5 13,7 14,8 15,1 14,3 15,2 15,2 13,8
14,75 12,35 15,15 14,85 14,35 14,85 15,50 13,80
Tabulka č. 5: Výsledky pokusů
16
Sestavením tabulky skončily přípravné a experimentální práce. Dále se budeme zabývat výpočty, jejichž cílem bude stanovit, které z faktorů ovlivňují významným způsobem dobu letu vrtulníku Y . Pro určení optimální úrovně faktorů a pro sestavení modelu dále potřebujeme vědět, které dvojice faktorů mají vzájemně významnou interakci, proto budeme také počítat jejich vliv na Y. V následující tabulce je uvedeno vzájemné působení faktorů, kde znaménka ve sloupcích SD, DU , SU , SDU získáme jako součin znamének v příslušných sloupcích.
pokus
S
D
U
SD
SU
DU
SDU
Y
1 2 3 4 5 6 7 8
− + − + − + − +
− − + + − − + +
− − − − + + + +
+ − − + + − − +
+ − + − − + − +
+ + − − − − + +
− + + − + − − +
14,75 12,35 15,15 14,85 14,35 14,85 15,50 13,80
Tabulka č. 6: Interakce faktorů
3.1.1. Výpočet efektu faktorů Následně musíme provést výpočet tzv. efektů (vlivů) jednotlivých faktorů. Efekt faktoru znamená přechod jednotlivých faktorů z dolní úrovně (-) na horní úroveň (+), který způsobí změnu doby letu Y . Pro výpočet efektu můžeme použít více metod. My se zaměříme na znaménkovou metodu a Yatesův algoritmus.
17
Znaménková metoda: Sečteme hodnoty ve sloupci Y, přičemž každá hodnota má takové znaménko, jaké znaménko má příslušný faktor v příslušném řádku, a tento součet vydělíme n/2, kde n =počet pokusů. Např.:
S=
−14,75 + 12,35 − 15,15 + 14,85 − 14,35 + 14,85 − 15,50 + 13,80 = −0,975. 4
Postupně vypočteme všechny efekty faktorů a jejich interakcí a zapíšeme do tabulky.
pokus
S
D
U
SD
SU
1 2 3 4 5 6 7 8
− + − + − + − +
− − + + − − + +
− − − − + + + +
+ − − + + − − +
+ − + − − + − +
DU SDU + + − − − − + +
− + + − + − − +
Y 14,75 12,35 15,15 14,85 14,35 14,85 15,50 13,80
efekt
faktor
14,45 průměr −0,975 S 0,75 D −0,025 SD 0,35 U 0,375 SU −0,7 DU −1,075 SDU
Tabulka č. 7: Efekt faktorů - znaménková metoda Druhou metodou výpočtu efektu je Yatesův algoritmus. Zde se postupuje tak, že si vytvoříme pomocnou tabulku se sloupci (1), (2), (3). Do sloupce (1) postupně zapisujeme součet 1. a 2. řádku, potom 3. a 4. řádku, 5. a 6. řádku, 7. a 8. řádku, potom následuje rozdíl 2. a 1. řádku, potom 4. a 3. řádku, 6. a 5. řádku, 8. a 7. řádku. Pro výpočet používáme hodnoty ze sloupce Y. Stejným způsobem vytvoříme sloupec (2), akorát použijeme hodnoty ze sloupce (1), a podobně sloupec (3), kde použijeme hodnoty ze sloupce (2). Počet sloupců vždy odpovídá počtu faktorů.
18
Hodnotu 1. řádku sloupce (3) vydělíme n =počet pokusů a tím získáme průměrnou dobu letu vrtulníku. Ostatní řádky vydělíme n/2 a tak dostaneme efekty faktorů a jejich interakcí. Posloupnost řádků v posledním sloupci je pevně daná, nedá se měnit. Z toho také vyplývá seřazení faktorů v posledním sloupci Tab. 6.
Y
(1)
(2)
(3)
dělitel
efekt
faktor
14,75 12,35 15,15 14,85 14,35 14,85 15,50 13,80
27,1 30,0 29,2 29,3 −2,4 −0,3 0,5 −1,7
57,1 58,5 −2,7 −1,2 2,9 0,1 2,1 −2,2
115,6 −3,9 3,0 −0,1 1,4 1,5 −2,8 −4,3
8 4 4 4 4 4 4 4
14,45 −0,975 0,75 −0,025 0,35 0,375 −0,7 −1,075
průměr S D SD U SU DU SDU
Tabulka č. 8: Efekt faktorů - Yatesův algoritmus
3.1.2. Test významnosti efektu K testování významnosti efektu potřebujeme nejprve znát rozptyl odhadu efektu. Ten se vypočte pomocí vzorce
s2e =
4σ 2 , n
kde σ 2 je rozptyl Y a n je počet pokusů (včetně opakování). σ 2 zjistíme pomocí veličiny s2 , která se v případě opakovaných pokusů vypočte jako
s2 =
v1 s21 + v2 s22 + · · · + vk s2k , v1 + v2 + · · · + vk
kde vi = ni − 1, ni =počet opakování i–tého pokusu, s2i je rozptyl i–tého pokusu. 19
Rozptyl s2i se vypočte pomocí vzorce n
s2i
1 X = (Yi − Y )2 . n − 1 i=1
Jelikož máme v našem příkladu pouze dvě opakování, lze s využitím vzorce pro aritmetický průměr Y =
Y1 −Y2 2
tento vztah zjednodušit
n
s2i
2
1 X 1 X = (Yi − Y )2 = (Yi − Y )2 = n − 1 i=1 2 − 1 i=1 = (Y1 −
Y1 − Y2 2 Y1 − Y2 2 ) + (Y2 − ) = 2 2
=(
2Y1 − Y1 − Y2 2 2Y2 − Y1 − Y2 2 ) +( ) = 2 2
=(
Y2 − Y1 2 Y1 − Y2 2 ) +( ) = 2 2
=
2(Y1 − Y2 )2 (Y1 − Y2 )2 = . 4 2
Nyní můžeme začít testovat. Nejprve stanovíme nulovou a alternativní hypotézu, podle které budeme testovat významnost efektu H0 : efekt faktoru je bezvýznamný, H1 : efekt faktoru je významný. K testování použijeme statistiku
t=
efekt . se
Pro stanovení kritického oboru zvolíme hladinu významnosti α = 0,05 a zjistíme kritickou hodnotu testovacího kritéria tn1 +n2 +...+nk −n (α), kde n1 , ..., nk jsou počty opakování pokusů a n je počet pokusů bez opakování. Konkrétně ni = 2 a n = 8.
20
Počítáme: 1) kritickou hodnotu tn1 +n2 +...+nk −n (α) = t16−8 (0,05) = 1,86, 2) rozptyl s2 =
0,125 + 3,645 + 0,245 + 0,125 + 0,005 + 0,245 + 0,18 + 0 = 0,571 8
a s2e =
4s2 4.0,571 = = 0,143, n 16 se = 0,378.
Pro efekty i, kde statistika t splní nerovnost |ti | > tn1 +n2 +...+nk −n (α), zamítáme nulovou hypotézu. To znamená, že jsou významné. Pro přehlednost uspořádáme všechny výpočty do tabulky.
pokus
Y1
Y2
Y1 − Y2
s2i
efekt
t
1 2 3 4 5 6 7 8
15,0 11,0 15,5 14,6 14,4 14,5 15,8 13,8
14,5 13,7 14,8 15,1 14,3 15,2 15,2 13,8
0,5 −2,7 0,7 −0,5 0,1 −0,7 0,6 0
0,125 3,645 0,245 0,125 0,005 0,245 0,180 0
S = −0,975 D = 0,75 SD = −0,025 U = 0,35 SU = 0,375 DU = −0,7 SDU = −1,075
−2,579 1,984 −0,066 0,926 0,992 −1,851 −2,844
Tabulka č. 9: Testování významnosti efektu Z tabulky vyplývá, že faktory S, D a interakce SDU v absolutní hodnotě převyšují kritickou hodnotu tn1 +n2 +...+nk −n (α). 21
3.2. Poloviční plán Budeme vycházet z Příkladu 1 a použijeme některé výpočty z úplného faktorového plánu. Určíme si hlavní faktory S, D a jejich kombinaci U = SD. Zavedeme faktor I obsahující jen +. V Tab. 9 je ukázán úplný faktorový plán se zařazením jednotkového faktoru I.
pokus
I
S
D
U
SD
SU
DU
SDU
1 2 3 4 5 6 7 8
+ + + + + + + +
− + − + − + − +
− − + + − − + +
− − − − + + + +
+ − − + + − − +
+ − + − − + − +
+ + − − − − + +
− + + − + − − +
Tabulka č. 10: Úplný faktorový experiment s jednotkovým faktorem Efekty faktorů v úplném plánu, viz Tab. 8: faktor efekt
S
D
U
SD
SU
DU
SDU
−0,975 0,75 0,35 −0,025 0,375 −0,7 −1,075
Nyní zjistíme efekty faktorů v každé polovině úplného plánu. 3.2.1. První polovina plánu Z Tab. 10 vybereme ty pokusy, které odpovídají definiční rovnici I = SDU a připojíme sloupec „úplný plánÿ obsahující čísla pokusů v úplném faktorovém plánu.
22
pokus
I
S
D
U
SD
SU
DU
1 2 3 4
+ + + +
+ − − +
− + − +
− − + +
− − + +
− + − +
+ − − +
SDU úplný plán + + + +
2 3 5 8
Tabulka č. 11: První polovina plánu Jak je vidět v tabulce, některé faktory nebo kombinace faktorů se rovnají (mají stejná znaménka ve sloupcích), tj. S = DU , D = SU , U = SD. Po úpravě každé rovnosti s využitím zavedených vlastností operací s faktory vždy dostaneme definiční rovnici I = SDU .
S = DU
·S
⇒
I = SDU,
D = SU
·D
⇒
I = SDU,
U = SD
·U
⇒
I = SDU.
Z toho vyplývá, že nám postačí uvažovat tabulku obsahující pouze faktory S, D, a jejich kombinaci U a efekty vypočteme jen pro tyto tři faktory, použijeme znaménkovou metodu.
pokus
S
D
U
úplný plán
Y1
Y2
Y
1 2 3 4
+ − − +
− + − +
− − + +
2 3 5 8
11,0 15,5 14,4 13,8
13,7 14,8 14,3 13,8
12,35 15,15 14,35 13,80
Tabulka č. 12: První polovina plánu s vybranými faktory Z důvodu zaměnitelnosti se efekty uvádějí pro součet faktorů a interakcí jako S + DU =
12,35 − 15,15 − 14,35 + 13,80 = −1,675, 2 23
D + SU =
−12,35 + 15,15 − 14,35 + 13,80 = 1,125, 2
U + SD =
−12,35 − 15,15 + 14,35 + 13,80 = 0,325, 2
což ale neznamená, že na každou připadá polovina. 3.2.2. Druhá polovina plánu Budeme postupovat stejně jako u první poloviny, akorát zvolíme definiční rovnici I = −SDU .
pokus
I
S
D
U
SD
SU
DU
1 2 3 4
+ + + +
− + + −
− + − +
− − + +
+ + − −
+ − + −
+ − − +
SDU úplný plán − − − −
1 4 6 7
Tabulka č. 13: Druhá polovina plánu Opět vidíme, že se objevují stejné zaměnitelné dvojice, takže můžeme vytvořit tabulku a spočítat efekty.
pokus
S
D
U
úplný plán
Y
1 2 3 4
− + + −
− + − +
− − + +
2 3 5 8
14,75 14,85 14,85 15,50
Tabulka č. 14: Druhá polovina plánu s vybranými faktory Tentokrát se efekty uvádějí pro rozdíl faktorů a interakcí jako S − DU =
−14,75 + 14,85 + 14,85 − 15,50 = −0,275, 2
D − SU =
−14,75 + 14,85 − 14,85 + 15,50 = 0,375, 2 24
U − SD =
−14,75 − 14,85 + 14,85 + 15,50 = 0,375. 2
3.2.3. Shrnutí V první a druhé polovině plánu vidíme, že faktory a interakce mají odlišné efekty od úplného faktorového plánu. Je to proto, že faktory S, D, U jsou „kontaminoványÿ efektem interakce. Musíme určit, kolik z celkového efektu připadá na jednotlivé interakce v součtu, což získáme vyřešením rovnic v následující tabulce. S + DU = −1,675 D + SU = 1,125 U + SD = 0,325 S − DU = −0,275 D − SU = 0,375 U − SD = 0,375 S = −0,975
D = 0,75
U = 0,35
DU = −0,7
SU = 0,375
SD = −0,025
Ukazuje se, že efekty jednotlivých faktorů a interakcí jsou stejné jako u úplného faktorového plánu, z čehož vyplývá, že při snížení počtu pokusů na polovinu nedochází ke ztrátě informací a tedy ani ke změně výsledků.
4. Regresní model Mějme deterministicky určené hodnoty délky křídla s1 , . . . , sn , d1 , . . . , dn , u1 , . . . , un , a měřené hodnoty y1 , . . . , yn . Označme Y = (y1 , . . . , yn )0 , (viz [2]). Uvažujme model Y = β0 + β1 S + β2 D + β3 U + β4 S 2 + β5 D2 + + β6 U 2 + β7 S D + β8 DU + β9 U S + kde je bílý šum, t.j. ∼ Nn (0, σ 2 I). Odhad regresní parametrů najdeme pomocí známého parametru ˆ = (X0 X)−1 X0 Y, kde β
1, 1, X= 1,
s1 , d1 , u1 , s21 , d21 , u21 , s1 d1 , s1 u1 , d1 u1 , s2 , d2 , u2 , s22 , d22 , u22 , s2 d2 , s2 u2 , d2 u2 , , ... 2 2 2 sn , d1 , u1 , s1 , d1 , u1 , s1 d1 , s1 u1 , dn un , 25
varianční matice tohoto odhadu je dána vztahem ˆ = σ 2 (X0 X)−1 . var(β)
4.1. Regresní model Box–Behnken design 4.1.1. Vzorový příklad V knize [1] je uveden příklad spouštění vrtulníku, kdy je použit plán nazývaný Box-Behnken design, který vypadá takto:
pokus
A
B
C
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0 0 1 0 −1 −1 1 1 −1 0 1 0 0 0 −1
−1 0 0 −1 −1 0 1 0 1 0 −1 1 1 0 0
1 0 1 −1 0 1 0 −1 0 0 0 −1 1 0 −1
Tabulka č. 15: Box-Behnken design v kódovaných proměnných
Faktor A = {6,5; 7,5; 8,5} je šířka papíru, B = {9; 10; 11} délka papíru a C = {−90; 0; 90} startovací úhel. V tomto konkrétním příkladu naměřili hodnoty, které jsou uvedeny v následující tabulce.
26
pokus
A
B
C
Y
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
7,5 7,5 8,5 7,5 6,5 6,5 8,5 8,5 6,5 7,5 8,5 7,5 7,5 7,5 6,5
9 10 10 9 9 10 11 10 11 10 9 11 11 10 10
90 0 90 −90 0 90 0 −90 0 0 0 −90 90 0 −90
3,3 3,1 2,2 4,2 1,1 5,3 1,3 1,8 3,9 6,1 3,3 0,8 2,2 6,2 2,1
Tabulka č. 16: Faktory s úrovněmi a naměřenými hodnotami Box a Behnken vytvořili příslušný regresní model s parametry β0 , ..., β9 . Tento model uvedeme s již vypočítanými parametry.
Yb = −295,14 + 30,34A + 38,12B + 0,000139C − 1,25A2 − 1,48B 2 − − 0,000127C 2 − 1,2AB − 0,00778AC + 0,00639BC. Pomocí programu Matlab jsme vytvořili 3D graf, který vypadá stejně jako v knize [1].
27
Obrázek č. 2: Graf doby letu vrtulníku
Obrázek č. 3: Vrstevnice
28
4.1.2. Vlastní příklad Podle vzorového příkladu z knihy [1] jsme provedli stejný pokus s jinými hodnotami faktorů, abychom mohli porovnat výsledky měření a výsledný graf. Opět vytvoříme tabulku s faktory S = {8; 9,5; 11} šířka papíru, D = {13; 14,5; 16} délka papíru a U = {−90; 0; 90} startovací úhel. Y vyjadřuje průměrnou dobu letu vrtulníku.
pokus
S
D
U
Y
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
9,5 9,5 11 9,5 8 8 11 11 8 9,5 11 9,5 9,5 9,5 8
13 14,5 14,5 13 13 14,5 16 14,5 16 14,5 13 16 16 14,5 14,5
90 0 90 −90 0 90 0 −90 0 0 0 −90 90 0 −90
13,4 15,0 13,5 14,7 14,8 15,4 14,5 14,5 15,5 15,2 12,9 14,9 14,1 13,8 15,6
Tabulka č. 17: Box–Behnken design
Do programu Matlab jsme zadali příslušnou regresní funkci, jíž aproximujeme naměřená data. Regresní model vypadá následovně
Yb = β0 + β1 S + β2 D + β3 U + β4 S 2 + β5 D2 + + β6 U 2 + β7 S D + β8 DU + β9 U S + β10 S DU.
29
Obrázek č. 4: Aproximace doby letu vrtulníku Vidíme, že graf vypadá úplně jinak než graf ze vzorového příkladu, proto se v další části pokusíme pomocí zvýšení počtu měření přiblížit původnímu grafu.
30
4.2. Dodatečná měření K rozšíření počtu pokusů použijeme 15 pokusů z Box-Behnken designu a 8 pokusů z úplného plánu, což máme dohromady 23 pokusů, a opět vytvoříme 3D graf a pro lepší názornost i graf vrstevnic.
Obrázek č. 5: Vrstevnice
Obrázek č. 6: Aproximace doby letu
31
U vrstevnic se ukazuje, že kdybychom použili vrtulníky z papíru o šířce menší než 8 a délce mezi 10 a 20, pak bychom mohli získat maximální možnou dobu letu při daných faktorech. Byly vybrány faktory S = {4; 6; 8; 9; 10; 11}, D = {10; 13; 14; 16; 18; 20} a U zůstává stejný. Výběr hodnot a poté naměřené hodnoty jsou uvedeny v tabulce.
pokus
S
D
U
Y1
Y2
Y
1 2 3 4 5 6 7 8 9 10 11 12 13
4 4 4 4 4 6 11 8 9 6 8 10 11
18 14 16 20 13 10 10 10 10 20 20 20 20
−90 −90 0 90 90 0 −90 0 90 0 −90 90 90
19,4 18,5 19,3 22,1 17,7 14,7 13,4 13,1 14,6 18,0 16,3 14,7 14,7
19,1 16,7 19,5 22,0 17,5 14,4 12,2 13,1 14,0 17,9 16,7 14,2 15,0
19,25 17,6 19,4 22,05 17,6 14,55 12,8 13,1 14,3 17,95 16,5 14,45 14,85
Tabulka č. 18: Dodatečná měření Můžeme si všimnout, že rozsah doby letu vrtulníku je mnohem větší než v předchozích pokusech. V některých případech je skoro dvojnásobná. Pro lepší představu, jak jsme postoupili v experimentování, opět vytvoříme graf vrstevnic a 3D graf, ve kterých použijeme hodnoty z tabulek 5 (úplný plán), 17 a 18.
32
Obrázek č. 7: Vrstevnice
Obrázek č. 8: Aproximace doby letu
33
Je vidět, že doba letu vrtulníku se zvýší, ale nezískali jsme stejný graf jako v knize [1]. Ve vzorovém příkladu použili rozměry v palcích a my jsme použili své vlastní rozměry v centimetrech. V následující tabulce přepočteme hodnoty, abychom mohli porovnat výsledky. palce cm
6,5
7,5
8,5
9
10
11
16,51 19,05 21,59 22,86 25,4 27,94
Z Tab. 16 vzorového příkladu vyplývá, že nejdelší doba letu (6,2 s) je u vrtulníku s faktory S = 7,5 palců = 19,05 cm, D = 10 palců = 25,4 cm, U = 0◦ a nejkratší doba (0,8 s) u vrtulníku s faktory S = 7,5 palců = 19,05 cm, D = 11 palců = 27,94 cm, U = −90◦ . Doba letu se zvýšila o 675%. My jsme dospěli k výsledku, že zvýšíme dobu letu, jestliže budeme snižovat faktor S a zvyšovat faktor D, tj. jestliže rozdíl mezi délkou a šířkou papíru bude co největší, přičemž musí platit D > S. Pokud jsme totiž měli papír o rozměrech S = 4 a D = 20, mohli jsme si při experimentování všimnout, že vrtulník letí kolmo dolů a pravidelně točí vrtulemi. Na rozdíl od papíru s rozměry S = 11 a D = 10, kde vrtulník letěl po spirálovité trajektorii a celý rotoval, až skoro do vodorovné polohy. Maximální doby letu (22,05 s) jsme dosáhli s faktory S = 4, D = 20, U = 90◦ a minimální doby letu (12,35 s) s faktory S = 11, D = 13, U = −90◦ , takže doba letu se zvýšila o 78,54%. Při porovnávání výsledků si nelze nevšimnout velice nápadného rozdílu v procentech zvýšení doby letu. Z našich pozorování vyplývá, že vrtulník z nápadně úzkého a dlouhého papíru letí déle než z papíru nápadně širokého a krátkého. V příkladu z knihy [1] použili vrtulníky z papíru o nevelkých rozdílech mezi šířkou a délkou a přesto zvýšily dobu letu o tolik procent. Zde se nabízí otázka, zda autor nepoužil k názornosti příkladu fiktivní data, která neodpovídají skutečnosti. 34
5. Závěr Cílem práce bylo především seznámení se základními principy při plánování experimentu a s různými typy faktorových plánů. Druhým cílem bylo jejich uplatnění na praktické úloze věnované „výrobě papírových vrtulníkůÿ. Sledovaným znakem kvality byla doba letu. Užitím faktorových plánů při výrobě papírových vrtulníků jsem získala vhodný regresní model popisující závislost sledovaného znaku jakosti — doby letu vrtulníku na výrobních faktorech. Použitelnost modelu jsem ověřila dalšími měřeními. Těmito faktory byly šířka a délka papíru, ze kterého byl vrtulník vyroben a dále úhel startu vrtulníku. Doba letu prvního vrtulníku byla 13,25 s. Po aplikaci designu experimentu a regresního modelu jsem získala hodnoty výrobních parametrů vrtulníku. Průměrná doba letu vrtulníku vyrobeného s vypočtenými parametry byla 22,05 s. Dosáhla jsem tedy zlepšení sledovaného znaku jakosti o 66,42 %. Na dosažených výsledcích je vidět, že design experimentu je mocný nástroj pro vylepšování jakosti výroby. Dosažené zlepšení u skutečných výrobků by mohlo znamenat obrovské finanční úspory. Při tvorbě práce jsem si vyzkoušela použití statistických metod na konkrétní aplikaci, získala nové vědomosti z designu experimentu a také si vylepšila programovací dovednosti v programu Matlab.
35
Literatura [1] ALLEN, T. T.: Introduction to Engineering statistics and Six Sigma. Springer–Verlag. London. 2006. [2] ANDĚL J.: Statistické metody. Matfyzpress. Praha. 1998 [3] JAROŠOVÁ, E.: Navrhování experimentů. VŠE. Praha 1998. [4] MAROŠ B., TRÁVNÍČEK T.: Plánování experimentu. 5th International Conference Aplimat. Bratislava. 2006 [5] MONTGOMERY D. C.: Design and Analysis of Experiments. Wiley. 2000 [6] TAGUCHI G., CHOWDHURY S., WU Y.: Taguchi’s Quality Engineering Handbook. Wiley. New Jersey. 2005 [7] TOŠENOVSKÝ, J. - NOSKIEVIČOVÁ, D.: Statistické metody pro zlepšování jakosti. Ostrava. Montanex. 2000. [8] Design of Experiments http://www.wikipedia.org/ [online 15. 3. 2010] [9] DOE - statisticky navržený experiment http://www.interquality.cz/INTERN%C3%8DKURZY/Procesyakvalita/ tabid/64/Default.aspx [online 20. 3. 2010] [10] Vliv neortogonality plánu experimentu na statistickou korektnost modelu http://dsp.vscht.cz/konference matlab/MATLAB08/prispevky/ 073 moravka.pdf [online 1. 4. 2010] [11] Výroba vrtulníku http://quest.nasa.gov/space/teachers/rockets/act11ws2.html [online 5. 4. 2010]
36
