1
6. PODSTAWY PLANOWANIA EKSPERYMENTU
6.1. Pojęcie i rola badań doświadczalnych
Przez eksperyment rozumiemy badanie jakiegoś zjawiska, polegające na celowym
wywoływaniu tego zjawiska lub jego zmian oraz obserwacji i pomiarach, umożliwiających
wnioskowanie o jego właściwościach. Eksperyment przeprowadza się na drodze badań
doświadczalnych mających na celu poznanie informacji o faktach, obiektach, zjawiskach
bądź procesach. Na podstawie zgromadzonych informacji określa się model zjawiska, który
stanowi jego reprezentację w postaci użytkowej wyrażającej istotne jego cechy. Istnieją
modele lingwistyczne, fizyczne, matematyczne. Najczęściej stosowany model
matematyczny jest opisem zjawiska za pomocą liczb, zmiennych, zbiorów, funkcji, relacji
itd. Znajomość modelu matematycznego umożliwia przewidywanie przebiegu zjawiska lub
zachowania obiektu w różnych warunkach. Podstawy teoretyczne badań doświadczalnych
ujmuje dziedzina wiedzy nazywana teorią eksperymentu. Obejmuje ona następujące
zagadnienia:
1. planowanie eksperymentów,
2. metodykę modelowania matematycznego,
3. technikę przeprowadzania pomiarów,
4. analizę wyników pomiarów.
Zanim powstała teoria eksperymentu badania doświadczalne dotyczyły prostych
obiektów o jednej wielkości wejściowej. Eksperyment polegał na przeprowadzeniu
pomiarów wielkości wyjściowej y dla L arbitralnie wybranych wartości wielkości
wejściowej x oraz wyznaczeniu funkcji aproksymującej y = f(x). Wraz ze wzrostem
złożoności analizowanych obiektów, spowodowanym zwiększeniem liczby wielkości
wejściowych, badania doświadczalne realizowano w oparciu o dwie metody: metodę badań
kompletnych i monoselekcyjnych.
W metodzie badań kompletnych doświadczenie odbywało się w następujący sposób:
a) dla każdej zmiennej x
k
, k = 1, 2, ..., s wybierano L wartości równomiernie
rozmieszczonych w przedziale [x
k min
, x
k max
],
b) dla każdej kombinacji wartości wielkości wejściowych wykonywano pojedynczy
pomiar,
c) na podstawie przeprowadzonych pomiarów wyznaczano funkcję aproksymującą f(
⋅
).
Graficzna interpretacja tej metody na przykładzie poszukiwania ekstremum obiektu o
dwóch zmiennych wejściowych przedstawiona została na rys. 6.1.
Metoda badań kompletnych była skuteczną metodą identyfikacji funkcji obiektu badań
dla jednej lub dwóch zmiennych wejściowych. Adoptowanie tej metody do obiektów o
większej liczbie zmiennych prowadziło do eksplozji kombinatorycznej wynikającej z liczby
koniecznych do wykonania pomiarów:
n = L
S
(6.1)
Przykładowo dla s = 10, L = 10, oraz przyjmując czas trwania jednego pomiaru równy
1s uzyskujemy całkowity czas badań kompletnych wynoszący ponad 317lat!
W celu ograniczenia liczby wykonywanych pomiarów opracowano metodę badań
monoselekcyjnych, w której zastosowano procedurę właściwą obiektowi badań o jednym
wejściu w odniesieniu do obiektu o wielu wejściach. W metodzie tej dla każdej zmiennej x
k
wybierano L wartości równomiernie rozmieszczonych w przedziale [x
k min
, x
k max
]. Następnie
2
dokonywano pojedynczego wyboru (monoselekcji) kolejnych wartości x
k
i badano wpływ
tej wielkości na wielkość wyjściową y.
x
2
x
1
y
= f(
x
1
, x
2
)
y
max
Rys. 6.1. Plan kompletny dla dwóch zmiennych wejściowych
Równocześnie przyjmowano, iż wartości pozostałych wielkości wejściowych są stałe:
x
q
= const, q = 1, 2, ..., s ; q ≠ k. W ten sposób całkowicie ignorowano współzależności
między wielkościami wejściowymi i zamiast funkcji wielu zmiennych y = f(x
1
, x
2
, ..., x
k
, ..., x
s
)
uzyskiwano jedynie zbiór wielu funkcji jednej zmiennej y = f
k
(x
k
) dla arbitralnych wartości
pozostałych zmiennych (rys. 6.2).
x
1
x
2
y
x
1 opt
y
y = f( x
1
, x
2
)
y = f
1
( x
1
)
x
2
= const
y = f
2
(x
2
)
x
1
= x
1 opt
= const
x
2
x
1
x
2 opt
y
max
Rys. 6.2. Badania monoselekcyjne dla dwóch zmiennych wejściowych.
Poszukiwane ekstremum funkcji nie zostaje wyznaczone
3
W porównaniu do metody badań kompletnych osiągnięto znaczną redukcję liczby
koniecznych do przeprowadzenia pomiarów:
(
)
∑
=
−
+
=
s
k
L
L
n
2
1
(6.2)
Dla s = 10, L = 10, oraz przyjmując czas trwania jednego pomiaru równy 1s całkowity
czas badań monoselekcyjnych uległ skróceniu do zaledwie 91s.
Ze względu na ograniczenia metody kompletnej i monoselekcyjnej przy analizie
złożonych obiektów zaistniała potrzeba opracowania nowych metod badawczych. Nastąpił
rozwój teorii eksperymentu. Powstały metody planowania badań doświadczalnych, które
umożliwiły zwiększenie ilości i jakości uzyskiwanej informacji naukowej. Zmniejszeniu
uległa liczba koniecznych do przeprowadzenia pomiarów, a więc zredukowano koszty i
czas trwania badań.
Badania kompletne i monoselekcyjne pozostawiały swobodę wyboru wartości wielkości
wejściowych (punktów pomiarowych) dla których realizowano eksperyment, natomiast
wyniki pomiarów analizowano matematycznie dopiero po przeprowadzeniu doświadczenia.
Wykorzystując zasady teorii eksperymentu ustala się wstępnie cel i metodę analizy
wyników pomiarów, natomiast punkty pomiarowe generowane są na podstawie
określonych procedur matematycznych. Pojedynczy punkt pomiarowy, będący
s–wymiarowym wektorem wartości wielkości wejściowych, nazywany jest układem planu
eksperymentu, natomiast zbiór wszystkich punktów pomiarowych stanowi plan
eksperymentu. W zależności od celu badań i stosowanych metod analizy uzyskuje się
rozmieszczenie punktów, które pozwala na:
- uwypuklenie poszukiwanych cech obiektu np. liniowości, współzależności zmiennych
wejściowych, niezależności wielkości wyjściowej od wielkości wejściowych,
- wyznaczenie ekstremum globalnego funkcji obiektu badań,
- zmniejszenie
nakładu obliczeniowego przy identyfikacji modelu obiektu.
Jak wspomnieliśmy celem badań doświadczalnych jest zwykle wyznaczenie modelu
obiektu badań. W ogólności rozróżnia się identyfikację strukturalną polegającą na ustalaniu
struktury modelu i wyznaczeniu wartości jego parametrów, oraz identyfikację
parametryczną polegającą na ustaleniu wartości parametrów modelu przy a priori danej
strukturze modelu. Wyznaczenie dokładnego modelu badanego obiektu jest bardzo trudne
ze względu na oddziaływanie na rzeczywisty obiekt wielu trudnych lub niemożliwych do
zmierzenia zakłóceń losowych. Z tego względu w teorii eksperymentu stosuje się
identyfikację parametryczną, w której przy założonej strukturze modelu i określonych
danych wejściowych wyznacza się takie wartości parametrów, które zminimalizują
niedokładność modelu. Przyjmuje się postać funkcji aproksymującej model obiektu
(najczęściej wielomian algebraiczny) i wyznacza się wartości współczynników tego
wielomianu korzystając z metody regresji.
Planowanie eksperymentu odbywa się zatem według następującego scenariusza:
a) charakterystyka obiektu badań polegająca na sformułowaniu zagadnienia
wymagającego rozwiązania na drodze doświadczalnej, ustaleniu wielkości
charakteryzujących obiekt badań (wielkości wejściowe, wyjściowe, stałe i zakłócające),
oraz przyjęciu liczby poziomów zmiennych wejściowych, czyli wybranie wartości, które
mogą przyjmować zmienne wejściowe,
b) ustalenie celu badań doświadczalnych, którym może być:
- identyfikacja modelu obiektu badań,
- optymalizacja empiryczna – wyznaczenie ekstremum globalnego modelu obiektu,
4
- badania eliminacyjne – eliminacja wielkości wejściowych o nieistotnym znaczeniu,
c) generacja lub wybór planu eksperymentu najlepiej dostosowanego do obiektu badań
i przyjętego celu badań doświadczalnych,
d) realizacja pomiarów w oparciu o wybrany plan doświadczenia,
e) analiza danych empirycznych zmierzająca do osiągnięcia założonego celu badań
doświadczalnych,
f) analiza merytoryczna rezultatów przeprowadzonych badań polegająca na logicznej
ocenie zjawisk związanych z badanym obiektem,
g) sformułowanie wniosków poznawczych, praktycznych i rozwojowych z
przeprowadzonych badań.
W dalszej części niniejszego rozdziału omówione zostały podstawowe zagadnienia z
zakresu planowania i analizy eksperymentu. Przedstawione wiadomości zaczerpnięto ze
specjalistycznej literatury poświęconej metodologii prowadzenia badań empirycznych [1-5],
adresowanej do inżynierów różnych specjalności, zajmujących się doświadczalnictwem
oraz identyfikacją modeli matematycznych obiektów.
6.2. Charakterystyka obiektu badań
Metody planowania eksperymentu mają charakter uniwersalny, niezależny od
merytorycznego obszaru, w którym realizowany jest eksperyment. Uzyskano to dzięki
wprowadzeniu uniwersalnego, przyczynowo-skutkowego modelu obiektu badań (rys. 6.3),
który opisują następujące wielkości modelujące:
a) niezależne wielkości wejściowe: {x
1
, x
2
, ..., x
k
, ..., x
s
}, k = 1, 2, ..., s,
b) zależne wielkości wyjściowe: {y
1
, y
2
, ..., y
p
, ..., y
w
}, p = 1, 2, ..., w,
c) wielkości stałe {c}, które mają wpływ
na działanie układu, ale ich wartości są
ustalone i niezmienne w czasie, przez
co mogą zostać pominięte w analizie
statystycznej,
d) wielkości zakłócające {h}, których
istnienie spowodowane jest oddziały-
waniem losowych czynników w
obiekcie badań na wielkości
wyjściowe oraz niedoskonałościami
metod i środków pomiarowych.
Do wielkości wejściowych wlicza się te, których wpływ na wielkości wyjściowe
interesuje realizatora badań. Mogą to być wielkości: fizyczne, chemiczne, techniczne,
ekonomiczne i inne. W odniesieniu do każdej wielkości wejściowej określa się zakresy
wartości x
k min
≤ x
k
≤ x
k max
, k = 1, 2, ..., s zakładając, iż eksperyment jest fizycznie
realizowalny jeżeli wartości mieszczą się w wyznaczonych przedziałach.
W celu uniezależnienia się od fizycznej interpretacji oraz zmniejszenia błędów
numerycznych podczas analizy danych empirycznych wielkości wejściowe normuje się do
bezwymiarowego przedziału [-
α, α]. Wielkość α jest nazywana ramieniem gwiezdnym i
stanowi jeden z elementów charakteryzujących plan eksperymentu.
Normowanie wielkości wejściowej x
k
realizowane jest zgodnie ze wzorem:
(
)
α
α
−
−
−
=
min
max
min
2
k
k
k
k
k
x
x
x
x
x(
(6.3)
OBIEKT
BADAŃ
x
1
x
2
:
x
k
:
x
s
y
1
y
2
:
y
p
:
y
w
Rys. 6.3. Uniwersalny, statyczny
model obiektu badań
5
Rzeczywistą wartość wielkości wejściowej uzyskuje się przez zastosowanie
przekształcenia odwrotnego zwanego denormowaniem:
(
)(
)
min
min
max
2
k
k
k
k
k
x
x
x
x
x
+
−
+
=
α
α
(
(6.4)
Jeżeli wszystkie wielkości modelujące są niezależne od czasu wówczas mamy do
czynienia z obiektem statycznym. W przeciwnym razie, jeżeli przynajmniej jeden z
parametrów jest funkcją czasu, wówczas operujemy pojęciem obiekt dynamiczny. Metody
identyfikacji obiektów dynamicznych są dużo bardziej złożone bowiem wymagają
rozwiązywania równań różniczkowych i w związku z tym nie będą omawiane w tym
podręczniku.
Analizę obiektu badań opisanego liczbą w > 1 wielkości wyjściowych sprowadza się
przez dekompozycję do analizy w obiektów o jednym wyjściu. Dalsze rozważania będą
dotyczyły obiektu zawierającego i-wejść x
k
, k = 1, 2, ..., s oraz jedno wyjście y.
6.3. Metody planowania eksperymentu
Kolejnym etapem badań doświadczalnych (po charakterystyce obiektu i ustaleniu celu
badań) jest wyznaczenie zbioru punktów pomiarowych czyli generacja lub wybór planu
eksperymentu. Algorytmy generacji planów ustala teoria eksperymentu na podstawie
określonych reguł matematycznych. Prowadząc badania doświadczalne najczęściej
wybieramy jeden plan z bazy planów eksperymentów dedykowany określonemu celowi
badawczemu np. identyfikacji modelu liniowego, optymalizacji modelu liniowo-
kwadratowego, badaniu wpływu poszczególnych składników na właściwości mieszaniny
chemicznej.
W zapisie matematycznym plan eksperymentu stanowi macierz
=
ns
nk
n
n
us
uk
u
u
s
k
s
k
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
...,
,
...,
,
,
...
...,
...,
...,
...,
...,
...,
,
...,
,
,
...
...,
...,
...,
...,
...,
...,
,
...,
,
,
...,
,
...,
,
,
2
1
2
1
2
2
22
21
1
1
12
11
X
(6.5)
gdzie: n – liczba układów planu eksperymentu; s – liczba zmiennych wejściowych.
Wiersz x
u
= [x
u1
, x
u2
, ..., x
uk
, ..., x
us
] macierzy X stanowi układ planu eksperymentu.
Istnieją różne klasyfikacje planów w zależności od struktury modelu identyfikowanego
obiektu, parametrów wielkości modelujących oraz celu badań doświadczalnych.
Liczba poziomów L zmiennych wejściowych x
k
determinuje plan L-poziomowy
(dwupoziomowy, trójpoziomowy, wielopoziomowy).
W zależności od proporcji liczby punktów pomiarowych n oraz liczby parametrów P
identyfikowanego modelu wyróżnia się:
a) plany nienasycone, w których n > P,
b) plany nasycone, w których n = P.
Rząd planu doświadczenia uzależniony jest natomiast od stopnia wielomianu
aproksymującego model obiektu. Aproksymacja wielomianem algebraicznym stopnia i
wymaga zastosowania planu i-tego rzędu.
6
W zależności od wpływu czasu na wartości zmiennych wejściowych wyróżniamy dwa
rodzaje planów:
a) statyczne – wszystkie wielkości modelujące są niezależne od czasu,
b) dynamiczne – przynajmniej jedna z wielkości modelujących jest funkcją czasu.
Najistotniejszym kryterium podziału jest przyjęty cel badań doświadczalnych.
a) Weryfikację istotności wpływu wielkości wejściowych na wielkość wyjściową
przeprowadza się w oparciu o plan randomizowany, który wprowadza element
losowości do zbioru punktów pomiarowych.
b) Identyfikację modelu obiektu badań najlepiej zrealizować w oparciu o plan
zdeterminowany, którego układy determinują ustalone założenia teoretyczne.
c) Przy wyznaczaniu ekstremum funkcji korzystamy z planu optymalizacyjnego.
Największe praktyczne zastosowanie mają plany zdeterminowane, a wśród nich plany
nazywane ułamkowymi lub poliselekcyjnymi.
Zdeterminowane metody planowania opisane zostały w dalszej części tego rozdziału.
Dokonano przy tym podziału w zależności od liczby poziomów zmiennych wejściowych.
Jako szczególny przypadek planów zdeterminowanych, mających zastosowanie w badaniu
zjawisk chemicznych, opisana została metoda planowania sympleksowego. Omówiono
także technikę planowania optymalnego.
6.3.1. Planowanie dwupoziomowe
Metoda planowania dwupoziomowego wykorzystywana jest do identyfikacji liniowych
modeli obiektów. W metodzie tej każda zmienna wejściowa przyjmuje tylko dwie wartości
(poziomy). Łączna liczba układów planu eksperymentu w planie dwupoziomowym wynosi
2
S
, gdzie s jest liczbą zmiennych wejściowych. Taki plan nosi nazwę dwupoziomowego
całkowitego lub kompletnego i oznaczany jest symbolem 2
S
. Dla większej liczby
zmiennych wejściowych przeprowadzenie eksperymentu całkowitego jest praktycznie
niemożliwe ze względu na dużą liczbę koniecznych do przeprowadzenia pomiarów.
W takim przypadku stosowane są plany ułamkowe (oznaczane symbolem 2
S-P
) zawierające
pewną liczbę układów z planu eksperymentu całkowitego. Istnieją zatem plany połówkowe,
ćwiartkowe, ósemkowe itd.
Unormowane zmienne wejściowe w planie dwupoziomowym przyjmują wartości
k
x(
={ -1, 1}. Rozmieszczenie układów planu całkowitego i ułamkowego (połówkowego)
dla trzech zmiennych wejściowych we współrzędnych unormowanych podano na rys. 6.4.
1
1
1
1
x(
2
x(
3
x(
1
1
1
1
x(
2
x(
3
x(
0
0
a) b)
Rys. 6.4. Rozmieszczenie we współrzędnych unormowanych układów planu dwupoziomowego:
a) całkowitego, b) ułamkowego (połówkowego)
7
Wybór konkretnych układów nie może być przypadkowy, gdyż plan ułamkowy dla
standaryzowanych zmiennych wejściowych powinien spełniać warunki:
a) symetrii układów względem środka eksperymentu,
b) ortogonalności, polegającej na zerowaniu wszystkich iloczynów skalarnych wektorów
kolumnowych unormowanej macierzy X,
c) równości sum kwadratów we wszystkich kolumnach unormowanej macierzy X.
6.3.2. Planowanie trójpoziomowe
Metoda planowania trójpoziomowego umożliwia identyfikację kwadratowego modelu
obiektu. W metodzie tej unormowane zmienne wejściowe przyjmują wartości
k
x(
={ -1, 0, 1}.
W planie trójpoziomowym całkowitym (oznaczenie 3
S
) występuje bardzo gwałtowny
wzrost liczby układów wraz ze wzrostem liczby wejść obiektu, stąd planowanie
trójpoziomowe ma bardzo małe możliwości praktyczne. Identyfikacja modelu na podstawie
planu trójpoziomowego jest bardziej skomplikowana niż na podstawie planu
dwupoziomowego. Z tego powodu nie zostały opracowane plany eksperymentów
trójpoziomowych ułamkowych.
6.3.3. Planowanie wielopoziomowe
Planowanie wielopoziomowe zapewnia identyfikację modeli liniowo-kwadratowych.
Szczególny przypadek planowania wielopoziomowego – planowanie pięciopoziomowe –
stanowi rozszerzenie planowania dwupoziomowego i jest najczęściej wykorzystywane w
praktyce.
Wyróżniamy trzy zasadnicze typy planowania wielopoziomowego:
1) planowanie kompozycyjne,
2) planowanie ortogonalne,
3) planowanie rotabilne.
Planowanie kompozycyjne jest rozwinięciem planowania dwupoziomowego typu 2
S
lub
2
S-P
o dwa rodzaje układów (dla zmiennych unormowanych):
a) gwiezdne typu (0, ..., 0,
±α, 0, ..., 0), w których zmieniane są kolejno wartości
zmiennych wejściowych między poziomami
±α dla pozostałych zmiennych na
poziomie 0, przy czym wielkość
α stanowi ramię gwiezdne planu,
b) centralne typu (0, 0, ..., 0) stanowiące centrum planu eksperymentu.
Przykłady planów kompozycyjnych dla unormowanych zmiennych wejściowych
przedstawia rys. 6.5.
1
1
1
1
x(
2
x(
3
x(
0
α
α
α
1
1
1
x(
2
x(
0
α
α
a)
b)
Rys. 6.5. Plany wielopoziomowe (kompozycyjne) we współrzędnych unormowanych dla różnej
liczby zmiennych wejściowych: a) dwóch, b) trzech
8
Liczba układów planu kompozycyjnego wynosi n = 2
S
+ 2s + 1 (2
S
układów planu
dwupoziomowego, 2s układów gwiezdnych i jeden układ centralny). Stąd podstawową
zaletą tego planu jest znaczne ograniczenie liczby układów w porównaniu do planu
trójpoziomowego, w szczególności dla większych wartości s.
Dobierając odpowiednią wartość ramienia gwiezdnego
α w planie kompozycyjnym
oraz zwiększając do n
o
liczbę układów w centrum planu można spełnić postulat
ortogonalności planu doświadczenia. Uzyskujemy w ten sposób znaczne uproszczenie
obliczeń przy wyznaczaniu parametrów modelu identyfikowanego obiektu oraz ocenie
statystycznej otrzymanych współczynników. Plan spełniający postulat ortogonalności
nazywany jest planem ortogonalnym.
Wartość ramienia gwiezdnego planu ortogonalnego dla określonych wartości s i n
o
wyznaczamy ze wzoru
(
)
[
]
d
d
d
n
n
S
n
n
−
+
+
=
0
2
2
1
α
(6.6)
gdzie: n
d
= 2
S
dla panu całkowitego lub n
d
= 2
S-P
dla planu ułamkowego.
Planowanie rotabilne ma na celu spełnienie postulatu niezależności planu od obrotu
układu współrzędnych w przestrzeni wielkości wejściowych. Zastosowanie planu
rotabilnego umożliwia identyfikację modelu o wariancjach zależnych tylko od odległości
od punktu centralnego eksperymentu.
Warunek rotabilności planu jest spełniony, jeżeli wartość ramienia gwiezdnego wynosi
4
2
S
=
α
dla panu całkowitego lub
4
2
P
S
−
=
α
dla panu ułamkowego.
W tablicy 6.1. podano zestawienie wartości ramienia gwiezdnego, zalecaną liczbę
układów w centrum planu n
o
oraz łączną liczbę układów n planu rotabilnego dla liczby
zmiennych wejściowych od 2 do 6.
Tablica 6.1.
Zestawienie optymalnych wartości parametrów dla planu rotabilnego
s
2 3 4 5 6
α
1,414 1,682 2,000 2,378 2,828
n
o
5 6 7 10 15
n
13 20 31 52 91
6.3.4. Planowanie sympleksowe
Specyficzny sposób planowania eksperymentu stosowany jest w przypadku badania
właściwości mieszaniny zależnej od jej składu.
Skład mieszaniny możemy opisać za pomocą wektora zmiennych x = [x
1
, x
2
, ..., x
s
],
przy czym na zmienne te są narzucone następujące ograniczenia fizykalne:
a) objętość mieszaniny jest stała, niezależna od jej składu:
∑
=
=
s
k
k
x
1
const (6.7)
b) zawartość każdego składnika w mieszaninie jest niezerowa:
x
k
≥ 0 , k = 1, 2, ..., s (6.8)
9
Przykładowo dla mieszaniny zawierającej trzy składniki wszystkie możliwe składy
mieszaniny leżą w przestrzeni trójwymiarowej na płaszczyźnie (warunek 6.7), ograniczonej
trójkątem (warunek 6.8). Sytuacja taka została zilustrowana na rys. 6.6.
x
1
x
2
x
3
a) b)
x
1
=1
x
2
=0
x
3
=0
x
1
=0
x
2
=1
x
3
=0
x
1
=0
x
2
=0
x
3
=1
x
1
=1/2
x
2
=1/4
x
3
=1/4
Rys. 6.6. Płaszczyzna składów mieszaniny o trzech składnikach:
a) w przestrzeni trójwymiarowej, b) w przestrzeni dwuwymiarowej
W ogólnym przypadku przestrzeń dopuszczalnych wartości s składników mieszaniny
jest sympleksem o s wierzchołkach na (s – 1)–wymiarowej hiperpłaszczyźnie.
W omawianych dotąd metodach planowania wartość wyjściową obiektu opisanego s
zmiennymi wejściowymi traktowano jako funkcję s zmiennych niezależnych
aproksymowaną za pomocą wielomianu algebraicznego stopnia R. W metodzie planowania
sympleksowego natomiast zmienne wejściowe związane są zależnością sumacyjną (6.7),
która zmniejsza liczę zmiennych niezależnych do s–1. Do opisu obiektu na który nałożone
są ograniczenia (6.7) i (6.8) stosuje się wielomian zredukowany stopnia R–1, który uzyskuje
się drogą odpowiednich przekształceń wielomianu algebraicznego stopnia R. Często
zachodzi potrzeba ograniczenia dużej liczby współczynników wielomianów
zredukowanych. Usuwając część z nich uzyskuje się uproszczone wielomiany
zredukowane.
Z uwagi na konieczność spełnienia warunków (6.7) i (6.8) sympleksowy plan
eksperymentu może zawierać wyłącznie układy stanowiące punkty leżące na sympleksie
s –wymiarowym.
Do wyznaczenia K współczynników wielomianu zredukowanego stopnia R dla s
zmiennych stosowane są plany siatkowe całkowite typu {s, R}. Plan siatkowy całkowity
{s, R} jest zbiorem układów określonych wzorami:
R
R
R
R
x
k
,...,
2
,
1
,
0
=
oraz
1
1
∑
=
=
s
k
k
x
(6.9)
Na rys. 6.7 podano przykładowe plany typu {3, R} mieszanin trójskładnikowych
s = 3 stopnia R = 1, 2, 3, 4.
x
1
=1
a) b)
c)
x
2
=1
x
3
=1
x
1
=1
x
2
=1
x
3
=1
x
1
=1
x
2
=1
x
3
=1
{3,1} {3,2}
{3,3}
Rys. 6.7. Całkowity plan sympleksowy dla trzech zmiennych stopnia:
a) pierwszego, b) drugiego, c) trzeciego
10
x
1
=1
x
2
=1
x
4
=1
{4,3}
x
3
=1
Rys. 6.8. Całkowity plan sympleksowy
trzeciego stopnia dla czterech zmiennych
Liczba układów planu całkowitego
siatkowego wynosi:
−
+
=
R
s
R
n
1
(6.10)
Rys. 6.8 przedstawia rozmieszczenie
układów planu typu {4, 3} mieszaniny
czteroskładnikowej.
Zasadniczą wadą planów siatkowych
całkowitych typu {s, R} jest badanie
właściwości mieszaniny przede wszystkim
na granicy sympleksu.
Wadę tę usuwa się stosując plany siatkowe ułamkowe typu {s}, stanowiące zbiór
układów określonych wzorami:
(1, 0, ..., 0), (
0
,...,
0
,
2
1
,
2
1
), (
s
s
s
s
1
,...,
1
,
1
,
1
) (6.11)
Przykład planu siatkowego ułamkowego przedstawia rys. 6.9.
x
1
=1
x
2
=1
x
4
=1
{4}
x
3
=1
x
1
=1
x
2
=1
x
3
=1
{3}
a) b)
Rys. 6.9. Ułamkowy plan sympleksowy dla trzech (a) i czterech zmiennych (b)
Liczba układów planu siatkowego ułamkowego wynosi:
1
2
−
=
s
n
(6.12)
Plany siatkowe typu {s} stosowane są dla liczby składników s ≥ 5 ze względu na
znacznie mniejszą liczbę współczynników wielomianów zredukowanych w odróżnieniu od
planów typu {s, R}.
6.3.5. Planowanie optymalne
Wraz z rozwojem technik obliczeniowych wymagania stawiane planom eksperymentu
polegające na ułatwianiu obliczeń związanych z wyznaczaniem współczynników modelu
stały się mało istotne. Postulowano natomiast opracowanie metod planowania
zapewniających niezależne szacowanie współczynników modelu, ogólniejsze stawianie
zadań planowania eksperymentu oraz zastosowanie do rozwiązywania tych zadań aparatu
analizy matematycznej w postaci teorii miary.
Nastąpił rozwój metod planowania optymalnego, w których analizując właściwości
macierzy informacyjnej Fishera M (patrz punkt 6.5.1 – Identyfikacja modelu obiektu
11
badań, wzór 6.33) optymalizowano wartości wariancji określonych parametrów w analizie
modelu metodą regresji. Najpopularniejsze plany optymalne to plany: D-, E-, A-, G- oraz V-
optymalne.
Plan D-optymalny minimalizuje wartość uogólnionej wariancji parametrów modelu.
Zadanie to sprowadza się do maksymalizacji wyznacznika macierzy informacyjnej M.
Plan E-optymalny minimalizuje wartość największej wariancji parametrów modelu,
poprzez maksymalizację najmniejszej wartości własnej macierzy M.
Plan A-optymalny minimalizuje średnią wariancję parametrów modelu. Zadanie to
polega na minimalizacji śladu macierzy kowariancyjnej C, będącej odwrotnością macierzy M.
Plany G-(V-) optymalne polegają natomiast na minimalizacji największej (średniej)
wartości wariancji prognozowanej wartości funkcji modelu (wzór 6.15) w punktach
stanowiących plan eksperymentu.
Z uwagi na szczególny charakter każdego przeprowadzanego eksperymentu wybór
konkretnego planu pozostaje zawsze w gestii badacza, który często kieruje się własnymi
kryteriami lub intuicją. Plan doświadczenia musi jednak spełniać podstawowe kryteria:
informatywności, realizowalności i efektywności.
Kryterium informatywności planu polega na jego zdolności do dostarczenia wymaganej
ilości informacji potrzebnej do osiągnięcia założonego celu badań.
Jeżeli np. celem badań jest wyznaczenie wielomianu aproksymującego stopnia i, o
liczbie współczynników niewiadomych P, wówczas powinny być spełnione zależności:
n ≥ P (6.13)
L ≥ i + 1 (6.14)
gdzie: n – liczba układów planu; L – liczba poziomów każdej wielkości wejściowej.
Zastosowanie tego kryterium do wyboru planu doświadczenia spowoduje odrzucenie
planów nie spełniających warunków (6.13) i (6.14).
Kryterium realizowalności planu polega na sprawdzeniu, czy analizowany plan jest
możliwy do fizycznej realizacji na stanowisku badawczym, w szczególności czy badany
obiekt będzie funkcjonował dla zadanego zbioru punktów pomiarowych. Plany nie
możliwe do realizacji są odrzucane, a o ostatecznym wyborze decyduje kryterium
efektywności.
Kryterium efektywności dotyczy ograniczenia kosztów i czasu badań poprzez
zmniejszenie liczby wykonywanych pomiarów. Ograniczenie liczby pomiarów zmniejsza
jednak możliwość dokładnego wyznaczenia funkcji obiektu badań. Trzeba zatem pójść na
kompromis między wzrostem dokładności identyfikacji funkcji obiektu uzyskiwanej przy
zwiększaniu liczby pomiarów, a redukcją kosztów badań uzyskiwaną z kolei przy
zmniejszaniu liczby pomiarów.
6.4. Realizacja pomiarów
Dla każdego układu x
u
zawartego w macierzy planu eksperymentu X wykonywany jest
pomiar, którego wynik stanowi u-element wektora kolumnowego wielkości wyjściowych
y = [y
1
, y
2
, ..., y
u
, ..., y
n
]
T
.
W rzeczywistym obiekcie badań na wielkość wyjściową y oprócz wielkości
wejściowych x
k
wpływają również wielkości zakłócające {h} mające charakter losowy.
Wielkość wyjściowa y jest zatem zmienną losową, do której opisu stosuje się dwie miary:
12
a) położenia (np. wartość oczekiwana),
b) rozproszenia (np. odchylenie standardowe).
Wyznaczenie wartości parametrów statystycznych zmiennej losowej wymaga
przeprowadzenia dodatkowych pomiarów (powtórzeń) w celu uzyskania próby z populacji
generalnej.
Przyjmuje się następujące warianty realizacji powtórzeń:
1. Dla każdego układu planu doświadczenia wykonuje się jednakową liczbę r powtórzeń.
Wariant ten stosuje się wówczas, gdy na podstawie analizy obiektu badań można
przypuszczać, że zakłócenia losowe zależą od wartości x
k
.
2. Jeżeli w planie doświadczenia występują jednakowe układy (np. plan ortogonalny,
rotalny) wówczas dla każdego układu wykonuje się jeden pomiar, a powtórzenia
uzyskuje się dzięki pomiarom wartości wielkości wyjściowej dla jednakowych
układów.
3. Powtórzenia realizuje się niezależnie od planu doświadczenia dla arbitralnie wybranych
wartości x
k
.
6.5. Analiza danych empirycznych
6.5.1. Identyfikacja modelu obiektu badań
Model obiektu badań przedstawia się w postaci zależności matematycznej nazywanej
funkcją modelu, opisującej związek między wielkościami wejściowymi x, a wielkością
wyjściową modelu
yˆ
:
yˆ
= f (x ; a) (6.15)
gdzie
a = [a
1
, a
2
, ...a
P
]
T
(6.16)
jest wektorem P parametrów.
Istotną kwestią jest ustalenie struktury modelu, a więc odpowiedni wybór funkcji
modelu. Przy złym wyborze nie uzyskamy dostatecznego dopasowania wyników pomiarów
do wartości wyjściowych modelu dla punktów pomiarowych nie wchodzących w skład
planu doświadczenia. Przyjęty model nie będzie wówczas adekwatny do obiektu badań i
uniemożliwi przewidywanie przebiegu zjawiska lub zachowania obiektu w różnych
warunkach.
Najczęściej stosowaną funkcją modelu jest liniowa kombinacja funkcji bazowych
yˆ
= a
0
f
0
(x) + a
1
f
1
(x) + … + a
P
f
P
(x), (6.17)
którą można zapisać w postaci wektorowej
yˆ
= [f(x)]a (6.18)
gdzie: f(x) = [f
0
(x), f
1
(x), ..., f
P
(x)] jest wektorem funkcji bazowych.
Jeżeli liczba pomiarów n wykonywanych w trakcie eksperymentu jest równa liczbie
parametrów P (plan nasycony) wówczas identyfikacja parametryczna modelu opisanego
funkcją (6.15) polega na rozwiązaniu układu P równań
f [x
1
(u), x
2
(u), ..., x
i
(u); a
1
, a
2
, ..., a
P
] = y(u), u = 1, 2, ..., n (6.19)
względem parametrów a
1
, a
2
, ..., a
P
.
13
Jeżeli liczba pomiarów n jest większa od liczby parametrów P (plan nienasycony)
wówczas identyfikację parametryczną przeprowadza się metodą regresji polegającą na
znalezieniu takich wartości parametrów a
1
, a
2
, ..., a
P
, dla których funkcja modelu
aproksymuje, najlepiej w sensie przyjętego kryterium, wyniki pomiarów y
przeprowadzonych dla ustalonego planu doświadczenia X.
Najbardziej rozpowszechnioną metodą aproksymacji jest metoda najmniejszych
kwadratów. W metodzie tej funkcjonałem podlegającym minimalizacji jest suma
kwadratów błędów aproksymacji
( )
[
]
( ) ( )
[
] (
) (
)
)
(
ˆ
ˆ
ˆ
T
2
2
a
y
y
y
y
ϕ
=
−
−
=
−
=
∆
=
∑
∑
u
y
u
y
u
S
(6.20)
gdzie:
y
ˆ jest wektorem kolumnowym wyjść modelu obliczonym ze wzoru (6.15), a
sumowanie wykonuje się względem numeru pomiaru u = 1, 2, ..., n.
Minimalizacja funkcjonału (6.20) sprowadza się do poszukiwania minimum funkcji
ϕ
zależnej od P zmiennych rzeczywistych a
1
, a
2
, ..., a
P
. Warunkiem koniecznym uzyskania
minimum funkcji
ϕ jest minimalizacja pochodnych cząstkowych
0
=
∂
∂
p
a
ϕ
; p = 1, 2, ..., P (6.21)
Warunek ten dla modelu liniowego względem funkcji bazy (6.17) jest również
wystarczający.
Dominującą postacią funkcji aproksymującej utworzonej na bazie modelu liniowego
względem funkcji bazy jest wielomian algebraiczny drugiego stopnia, z podwójnymi
iloczynami stanowiącymi tzw. interakcje, o ogólnej postaci:
...
...
...
ˆ
2
1
12
2
2
1
11
1
1
0
+
+
+
+
+
+
+
+
=
x
x
a
x
a
x
a
x
a
x
a
a
y
i
ii
i
i
(6.22)
Przykładowo dla dwóch zmiennych wejściowych x
1
i x
2
wielomian drugiego stopnia
przedstawia się następująco:
2
2
22
2
1
11
2
1
12
2
2
1
1
0
ˆ
x
a
x
a
x
x
a
x
a
x
a
a
y
+
+
+
+
+
=
(6.23)
Wielomian ten można zapisać w postaci wektorowej (6.18) przyjmując następujące
postacie funkcji bazowych:
f
1
(x) = 1 ; f
2
(x) = x
1
; f
3
(x) = x
2
; f
4
(x) = x
1
x
2
; f
5
(x) = x
1
2
; f
6
(x) = x
2
2
(6.24)
oraz zapisując wektor funkcji bazowych oraz wektor parametrów w postaci:
f(x) = [1 x
1
x
2
x
1
x
2
x
1
2
x
2
2
] , a = [a
0
a
1
a
2
a
12
a
11
a
22
]
T
(6.25)
W oparciu o macierz planu doświadczenia X oraz na podstawie zdefiniowanych funkcji
bazowych (6.24) tworzona jest macierz F o wymiarze n
× N
b
, która odgrywa rolę macierzy
obserwacji nowych zmiennych wejściowych. Liczba tych zmiennych jest równa liczbie
identyfikowanych parametrów N
b
. Tworzenie macierzy F polega na wykonaniu na
argumentach x
ku
macierzy X operacji określonych przez poszczególne funkcje bazowe. Dla
wielomianu drugiego stopnia macierz F ma postać:
14
=
2
2
2
1
2
1
2
1
2
2
2
1
2
1
2
1
2
22
2
21
22
21
22
21
2
12
2
11
12
11
12
11
,
,
,
,
,
1
...
...,
...,
...,
...,
,
1
,
,
,
,
,
1
...
...,
...,
...,
...,
,
1
,
,
,
,
,
1
,
,
,
,
,
1
n
n
n
n
n
n
u
u
u
u
u
u
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
F
(6.26)
Identyfikacja parametryczna modelu dla planu nasyconego sprowadza się do
rozwiązania układu równań
Fa = y (6.27)
Dla planu nienasyconego układ ten jest sprzeczny ze względu na błędy pomiarów
wielkości wyjściowej y i/lub nieodpowiedni dobór struktury modelu. W takiej sytuacji
wykorzystując warunek minimalizacji funkcjonału danego wzorem (6.20) rozwiązuje się
tzw. równania normalne określone wzorem
F
T
Fa = F
T
y (6.28)
Jeżeli w trakcie wykonywania pomiarów przyjęto liczbę powtórzeń u-tego układu planu
doświadczenia równą r
u
wówczas równania normalne zapisuje się w postaci
F
T
RFa = F
T
R
y
(6.29)
gdzie R jest macierzą diagonalną powtórzeń pomiarów o wymiarze n
× n:
R = diag (r
1
, r
2
, ..., r
n
) (6.30)
Wektor
y
jest natomiast wektorem wartości średnich wielkości wyjściowej dla
kolejnych układów planu doświadczenia
[
]
T
2
1
,...,
,...,
,
n
u
y
y
y
y
=
y
(6.31)
∑
=
=
u
r
j
j
u
u
u
y
r
y
1
)
(
1
(6.32)
Rozwiązanie układu równań (6.28) wymaga w pierwszej kolejności wyznaczenia
macierzy informacyjnej Fishera o postaci:
M = F
T
RF (6.33)
Następnie oblicza się wyznacznik macierzy informacyjnej i sprawdza warunek
det (M) ≠ 0 (6.34)
Jeżeli warunek ten nie jest spełniony wówczas macierz M jest osobliwa i należy
zmodyfikować plan doświadczenia. W przeciwnym razie oblicza się macierz odwrotną
(kowariancyjną) C = M
-1
i rozwiązuje układ równań normalnych ze wzoru
a = CF
T
R
y
(6.35)
15
Metoda najmniejszych kwadratów rozwiązywania układu równań normalnych jest
szczególnie przydatna do aproksymacji wyników pomiarów wielomianami liniowymi oraz
wielomianami drugiego stopnia. Dla wielomianów wyższych rzędów mogą wystąpić
trudności w ocenie zgodności funkcji modelu z zachowaniem się obiektu badań w pełnym
zakresie wartości wejściowych.
Do identyfikacji modeli obiektów, oprócz wielomianów algebraicznych, można
stosować również inne postacie funkcji aproksymujących:
- funkcje
eksponencjalne, np.
x
a
e
a
y
1
0
=
- funkcje
potęgowe, np.
1
0
a
x
a
y
=
Warunkiem zastosowania metody najmniejszych kwadratów jest wówczas
przeprowadzenie „linearyzacji” funkcji poprzez zastosowanie nowych zmiennych
pomocniczych.
Zaleca się przeprowadzanie aproksymacji dla unormowanych wartości wielkości
wejściowych. Zmniejszeniu ulegają wówczas błędy obliczeń numerycznych spowodowane
zaokrągleniami przy wielokrotnym wykonywaniu operacji matematycznych. Możliwa jest
również weryfikacja istotności współczynników funkcji testem t-Studenta.
Przedstawiona metoda rozwiązywania układu równań normalnych wymaga obliczenia
macierzy odwrotnej. Operacja odwracania macierzy jest jednak złożona obliczeniowo i
podatna na błędy numeryczne wynikające z zaokrągleń pośrednich wyników obliczeń. Z
tego względu opracowano szereg metod rozwiązywania układów równań nie
wymagających przeprowadzania operacji odwracania macierzy np. metoda eliminacji
Gaussa-Jordana, rozkład względem wartości szczególnych SVD, metody dekompozycji
LU, Choleskiego i QR. Każda z wymienionych metod odznacza się pewnymi
właściwościami, które predestynują ją do rozwiązania układu równań normalnych (6.28),
utworzonego na podstawie wybranego planu doświadczenia i z którym związana jest
określona struktura macierzy informacyjnej M.
Identyfikację modelu można zrealizować poprzez interpolację wyników pomiarów
funkcjami sklejanymi – splajnami (ang. spline). Metoda ta polega na interpolowaniu
kolejnych par danych {x
u
, y
u
} wielomianami algebraicznymi stopnia m w taki sposób, aby
uzyskać funkcję gładką tzn. ciągłą i mającą ciągłe pochodne do rzędu m – 1 włącznie na
całym przedziale [x
k min
, x
k max
]. Konieczne jest również spełnienie twierdzenia o
jednoznaczności funkcji sklejanej, orzekające że m – 1 pochodne dla argumentów x
k min
i x
k max
powinny wynosić zero. Zastosowanie metody interpolacji funkcjami sklejanymi jest
uzasadnione dla małych wartości błędów pomiarów wielkości wyjściowej y.
6.5.2. Weryfikacja adekwatności modelu obiektu
Model obiektu opisuje jego właściwości i zachowanie tylko w przybliżeniu.
Spowodowane jest to niedokładnością wyznaczenia parametrów modelu oraz
nieadekwatnością struktury modelu.
Na niedokładność wyznaczenia parametrów modelu mają wpływ następujące czynniki:
- błędy przyjętej metody identyfikacji parametrów modelu,
- błędy obliczeń numerycznych,
- błędy danych użytych do identyfikacji parametrów modelu.
Nieadekwatność struktury modelu wynika natomiast z trzech czynników:
- pominięcia wśród wielkości modelujących obiekt, czynników istotnych dla przebiegu
zjawisk w obiekcie,
- niewłaściwą specyfikacją wielkości modelujących obiekt,
16
- przyjęciem niewłaściwego typu równania modelu.
Oceny adekwatności modelu dokonuje się na dwa sposoby.
Pierwsza metoda polega obliczeniu wartości błędu aproksymacji wybraną funkcją f(
⋅)
i porównaniu jej z pewną arbitralnie wybraną wartością dopuszczalną e
d
. Jeśli obliczona
wartość błędu e
max
jest mniejsza od e
d
wówczas uznaje się wyznaczony model za
adekwatny.
Najczęściej stosuje się następujące definicje błędów aproksymacji:
a) maksymalny bezwzględny błąd aproksymacji
u
u
y
y
ˆ
max
max
−
=
ε
(6.36)
b) maksymalny błąd względny
%
100
ˆ
max
max
⋅
−
=
u
u
u
y
y
y
δ
(6.37)
c) błąd średniokwadratowy
(
)
n
y
y
n
u
u
u
∑
=
−
=
1
max
ˆ
δ
(6.38)
gdzie:
y
u
– wartość wielkości wyjściowej w u-układzie planu doświadczenia,
u
yˆ
– wartość wielkości wyjściowej obliczonej z wyznaczonej funkcji aproksymującej,
n – liczba układów planu doświadczenia.
Innym sposobem jest zastosowanie statystycznego testu istotności – testu F
(Snedecora), w którym weryfikuje się statystycznie hipotezę:
H:
σ
a
2
⇔ σ
2
, (6.39)
porównując wariancję błędów aproksymacji (wariancję adekwatności)
σ
a
2
z wariancją
niedokładności pomiarów wielkości wyjściowej
σ
2
.
Przyjmuje się następujące hipotezy:
1) hipoteza zerowa H
0
:
σ
a
2
=
σ
2
oznaczająca, iż model jest adekwatny,
2) hipoteza alternatywna: H
1
:
σ
a
2
>
σ
2
oznaczająca, iż model nie jest adekwatny.
Procedura weryfikacji statystycznej dla jednakowej liczby powtórzeń r we wszystkich
układach planu eksperymentu jest następująca:
a) Oblicza się wartość funkcji testowej
2
2
σ
σ
a
F
=
, (6.40)
przy czym:
(
)
∑
=
−
=
n
u
u
u
a
y
y
f
r
1
2
2
2
ˆ
σ
(6.41)
17
∑
=
−
=
n
u
u
f
r
1
2
1
2
1
σ
σ
,
(
)
1
1
2
)
(
2
−
−
=
∑
=
r
y
y
r
j
u
j
u
u
σ
(6.42)
W powyższych wzorach f
1
= n(r–1) jest liczbą stopni swobody wariancji
σ
2
, natomiast
f
2
= n – N
b
jest liczbą stopni swobody wariancji
σ
a
2
.
b) Na podstawie rozkładu F (Snedecora) odczytuje się z tablicy statystycznej wartość
krytyczną
1
2
,
,
f
f
F
α
odpowiadającą założonemu poziomowi ufności
α.
c) Sprawdza się warunek F ≤
1
2
,
,
f
f
F
α
. Jeśli warunek jest spełniony wówczas nie ma
podstaw do odrzucenia hipotezy zerowej i przyjmuje się, że model jest adekwatny. W
przeciwnym razie prawdziwa jest hipoteza alternatywna, czyli model nie jest
adekwatny.
Stwierdzenie na podstawie jednej z wymienionych metod nieadekwatności modelu
obiektu oznacza konieczność ponownego przeprowadzenia aproksymacji funkcją o innej
postaci lub zwiększenie liczby pomiarów dla każdego układu planu doświadczenia.
6.5.3. Weryfikacja istotności współczynników funkcji aproksymującej
Charakteryzując obiekt badań przyjmuje się określoną liczbę zmiennych wejściowych.
Nie ma jednak pewności czy wszystkie zdefiniowane zmienne wejściowe mają wpływ na
działanie obiektu. Stwierdzenie braku skorelowania określonej zmiennej wejściowej x
k
ze
zmienną wyjściową y umożliwia uproszczenie modelu badań poprzez usunięcie zmiennej
x
k
. Działanie takie jest uzasadnione głównie ze względów ekonomicznych, gdyż prostszy
model oznacza mniejszą ilość sprzętu technicznego nie-zbędnego do przeprowadzenia
pomiarów oraz uproszczenie obliczeń matematycznych.
Informacja o wpływie kolejnych wielkości wejściowych x
k
na wielkość wyjściową y jest
ukryta w wartościach współczynników funkcji aproksymującej. Przykładowo jeśli
wszystkie współczynniki przy x
2
wynoszą zero tzn. a
2
= a
22
= a
12
= 0 wówczas można
stwierdzić, że wielkość wyjściowa y nie zależy od wielkości wejściowej x
2
. Gdyby
natomiast współczynniki przy x
2
wynosiły: a
22
= a
12
= 0 oraz a
2
≠ 0 wówczas można
wyciągnąć wniosek, że wielkość x
2
wpływa na wielkość wyjściową, ale tylko liniowo. Jak
widać analiza współczynników funkcji aproksymującej jest bardzo istotna dla realizatora
badań, który uzyskuje w ten sposób istotne informacje o sposobie działania obiektu.
Analiza ta nosi nazwę weryfikacji istotności współczynników funkcji aproksymującej.
Realizowana jest w oparciu o test t-Studenta oraz ocenę wartości kowariancji wszystkich
par współczynników {a
i
, a
j
} funkcji aproksymującej f(
⋅). Kowariancję współczynników a
i
oraz a
j
oblicza się ze wzoru:
cov(a
i
, a
j
) = c
ij
⋅
σ
2
(6.43)
gdzie c
ij
jest elementem macierzy kowariancyjnej C stanowiącej odwrotność macierzy
informacyjnej M (wzór 6.31). Wykrycie nieistotnych współczynników funkcji
aproksymującej na podstawie testu t-Studenta lub ich wzajemnego skorelowania
(niezerowej wartości kowariancji) wskazuje na konieczność uproszczenia modelu. Po
wyznaczeniu funkcji aproksymującej należy ponownie przeprowadzić weryfikację jej
adekwatności. Dopiero pozytywne przejście tej weryfikacji jest podstawą eliminacji
nieistotnych współczynników.
18
6.5.4. Optymalizacja funkcji aproksymującej
Optymalizacja wyznaczonej funkcji aproksymującej f(
⋅
) polega na znalezieniu
ekstremum globalnego (maksimum lub minimum) w zadanych przedziałach wartości
wielkości wejściowych [x
k min
, x
k max
]. Istnieje wiele metod optymalizacji dedykowanych
określonym rodzajom funkcji, jednak żadna z metod nie ma charakteru uniwersalnego
gwarantującego uzyskanie ekstremum globalnego dowolnej funkcji wielu zmiennych.
Metody gradientowe optymalizacji bazują na obliczeniach wartości pochodnych
(cząstkowych, kierunkowych) optymalizowanej funkcji i wymagają spełnienia warunku
różniczkowalności funkcji w całym zakresie zmienności wartości wielkości wejściowych,
co przy pewnych postaciach funkcji aproksymujących nie jest spełnione. Ponadto
algorytmy te cechuje skłonność do utykania w ekstremach lokalnych bez gwarancji
osiągnięcia ekstremum globalnego. Z tego względu opracowano szereg algorytmów
bezgradientowych (np. metoda Monte Carlo, metoda symulowanego wyżarzania, metoda
wg Gaussa-Seidela), które w obliczeniach nie wykorzystują pochodnych optymalizowanej
funkcji. Algorytmy te są mniej efektywne od algorytmów gradientowych pod względem
dokładności wyznaczania ekstremum i cechuje je także znaczne wydłużenie procesu
obliczeniowego. W przeciwieństwie do metod gradientowych są jednak skuteczniejsze w
wyznaczaniu ekstremów globalnych.
Z uwagi na wady każdej z dostępnych metod optymalizacji zaleca się równoczesne
wykonywanie obliczeń z zastosowaniem dwóch lub trzech algorytmów oraz wybór
najlepszego wyniku.
6.6. Analiza merytoryczna wyników badań doświadczalnych
Istotnym krokiem realizacji eksperymentu jest analiza merytoryczna wyników
przeprowadzonych badań. Sposób jej realizacji jest uzależniony od natury obiektu badań.
Obiekt fizyczny wymaga na przykład odrębnej analizy niż obiekt chemiczny, a jeszcze
innej niż obiekt ekonomiczny. Niezależnie jednak od natury obiektu należy poddać ocenie
poprawność realizacji badań i skuteczność zastosowanych metod analizy danych
empirycznych.
Wyznaczenie funkcji aproksymującej modelu obiektu oraz przeprowadzenie jej
optymalizacji powinno być uzupełnione sporządzeniem odpowiednich wykresów, które
pozwolą zweryfikować stopień osiągnięcia założonych na wstępie celów eksperymentu
naukowego. Wykres funkcji aproksymującej można przedstawić na płaszczyźnie lub w
przestrzeni tylko dla jednej lub dwóch zmiennych wejściowych. Dla większej liczby
zmiennych sporządza się wykresy zależności zmiennej wyjściowej od każdej zmiennej
wejściowej z osobna, przy ustalonych wartościach pozostałych wielkości wejściowych.
Wartości te można wybrać arbitralnie lub przyjąć wartości optymalne, dla których funkcja
aproksymująca osiąga ekstremum. Uzyskujemy wówczas tzw. cięcia przestrzeni
wielowymiarowej w ekstremum globalnym.
Ostatecznie można sformułować wnioski poznawcze, praktyczne i rozwojowe z
przeprowadzonych badań, które mogą dotyczyć poznania nowych zjawisk, wdrożenia
nowych technologii oraz wskazać zagadnienia wymagające dalszego rozpoznania na drodze
doświadczalnej.
19
6.7. Inteligentne systemy planowania eksperymentu
Kluczową rolę w rozwoju nowoczesnych metod planowania eksperymentu zaczynają
odgrywać systemy inteligentne wywodzące się z obszaru inżynierii systemów
inteligentnych (ang. Intelligent Systems Engineering). Powstała koncepcja drugiej
generacji planów eksperymentu (2GD), tzw. planów inteligentnych eksperymentu (akronim:
InDE). Konieczność wprowadzania nowych metod planowania wynika z istotnych
ograniczeń dotychczasowych metod, określanych jako pierwsza generacja planów
eksperymentu (1GD).
Metodologia teorii eksperymentu w pierwszej generacji planowania eksperymentów
wymagała dopasowania przeprowadzanego eksperymentu do jednego z istniejących
planów. Plany drugiej generacji są natomiast generowane poprzez inteligentne
dostosowywanie do warunków realizacji eksperymentu.
Na gruncie nowych metod planowania dominujące dotychczas procedury weryfikacji
postawionych hipotez są wypierane przez metody eksploracyjne (Data Mining). Wśród
nich można wymienić metodę rekurencyjnego podziału (analizę dyskryminacyjną), metody
wizualne, sztuczne sieci neuronowe i algorytmy genetyczne.
Metody eksploracyjne oparte na sztucznych sieciach neuronowych stosowane do
aproksymacji funkcji badanego obiektu marginalizują aproksymację numeryczną. Metody
numeryczne wymagają bowiem dokonania wielu założeń, które często nie mają fizycznego
odzwierciedlenia w rzeczywistych obiektach. Można tu wymienić następujące
ograniczenia:
1. Konieczność realizacji eksperymentu dla ściśle określonego zbioru wartości
wejściowych występujących w wybranym planie eksperymentu, pomimo iż wartości te
mogą znacznie różnić się od obowiązujących standardów lub też nie być realizowalne w
praktyce,
2. Zachowanie stałej dla wszystkich wielkości wejściowych liczby wartości wielkości
wejściowych, co może odbiegać od fizycznych możliwości generacji poszczególnych
sygnałów wejściowych obiektu,
3. Rozmieszczenie układów planu eksperymentu zwykle z zachowaniem określonej
symetrii tworzącej regularny obszar wielkości wejściowych, co nie jest skutecznym
mechanizmem identyfikacji funkcji rzeczywistych obiektów.
Warunki te można odrzucić, jeżeli wyniki pomiarów uzyskane w oparciu o inteligentny
plan eksperymentu będą stanowiły dane do aproksymacji neuronowej. W ten sposób można
indywidualnie generować plan eksperymentu do założonych z góry cech właściwych
konkretnemu obiektowi badań i warunków realizacji eksperymentu. Istotna zaleta planu
inteligentnego polega na możliwości jego wielokrotnego generowania, aż do uzyskania
rezultatu spełniającego postulaty badacza.
6.8. Komputerowe wspomaganie badań doświadczalnych
Jeszcze do niedawna złożoność obliczeniowa stosowanych algorytmów
matematycznych stanowiła istotną przeszkodę w rozwoju teorii eksperymentu. Żmudne i
czasochłonne obliczenia wymuszały opracowywanie specjalnych planów eksperymentu
ułatwiających obliczenia oraz ograniczenie liczby przeprowadzanych analiz.
Obecnie dysponujemy znacznymi mocami obliczeniowymi jakie oferują nam
komputery osobiste i problem złożoności obliczeniowej algorytmów nie ma już istotnego
znaczenia. Szeroko rozpowszechnione jest oprogramowanie CADEX/DOE (akronim:
Computer Aided Design and Analysis of Experiments / Design of Experiments)
wspomagające prowadzenie badań doświadczalnych (tablica 6.2).
20
Tablica 6.2.
Najbardziej znane o
programowanie CADEX/DOE
nazwa programu
firma lub uczelnia oferująca program
Design Expert
Stat-Ease
DAX-Expert – System ekspertowy
planowanie i analiza eksperymentu
Politechnika Krakowska
STATISTICA – moduł Planowanie
Doświadczeń
StatSoft
STATGRAPHICS PLUS
Manugistics
Matlab – Statistics Toolbox
MathWorks
JMP – The Statistical Discovery Software
SAS Institute
Oferowane oprogramowanie CADEX/DOE zawiera obszerną pomoc elektroniczną
wyjaśniającą wiele spraw związanych z planowaniem eksperymentu i podającą wskazówki
ułatwiające realizację badań doświadczalnych. Obsługa programów wymaga jednak od
użytkownika gruntownego przygotowania merytorycznego oraz umiejętności interpretacji
prezentowanych wyników analiz.
Przykład 1
Rozważmy eksperyment polegający na identyfikacji modelu zjawiska ruchu pocisku
balistycznego wystrzelonego z wyrzutni z prędkością początkową v
0
i ustawionej pod
kątem
ϕ do powierzchni Ziemi. Analizujemy wpływ dwóch wielkości wejściowych:
x
1
= v
0
i x
2
=
ϕ na zasięg pocisku y = s, tj. dystans pomiędzy wyrzutnią a miejscem
zetknięcia się pocisku z Ziemią.
Wprowadzamy następujące ograniczenia na wartości zmiennych wejściowych:
1) zakres wysokości ustawienia działa: v
0min
= 100 m/s, v
0max
= 200 m/s,
2) zakres kąta ustawienia działa:
ϕ
min
= 0,175rad = 10
°,
ϕ
max
= 1,396rad = 80
°.
Celem badań jest wyznaczenie funkcji modelu procesu w postaci wielomianu
algebraicznego drugiego stopnia o postaci:
5
2
2
1
3
3
2
2
1
2
2
2
1
1
ˆ
x
x
a
x
x
a
x
x
a
y
+
+
=
(6.44)
oraz przeprowadzenie optymalizacji funkcji modelu polegającej na znalezieniu ekstremum
globalnego w zadanych przedziałach wartości. Na podstawie charakterystyki obiektu badań
i założonego celu badań wybrano kompozycyjny, centralny plan doświadczenia o wartości
ramienia gwiezdnego
α = 1,414 zawierający n = 9 układów (tablica 6.3). Dla każdego
układu wykonano r = 10 pomiarów zasięgu pocisku i wyznaczono wartości parametrów
statystycznych: wartość oczekiwaną
u
y
i wariancję
σ
u
2
(tablica 6.3).
Identyfikację funkcji modelu przeprowadzono rozwiązując układ równań normalnych
(6.35). Uzyskano w ten sposób następujące wartości parametrów:
a
1
= 0,203 a
2
= -0,130 a
3
= -0,020
Adekwatność przyjętej struktury modelu weryfikujemy za pomocą testu F – Snedecora,
przyjmując hipotezy:
1) zerową H
0
:
σ
a
2
=
σ
2
oznaczająca, iż model jest adekwatny,
2) alternatywną: H
1
:
σ
a
2
>
σ
2
oznaczająca, iż model nie jest adekwatny.
21
Tablica 6.3.
unormowane
wartości
zmiennych
wejściowych
rzeczywiste
wartości
zmiennych
wejściowych
wielkość
wyjściowa
wartość
obliczona
z funkcji
modelu
numer
układu
u
1
u
x(
2
u
x(
x
u1
x
u2
u
y
2
u
σ
yˆ
1 -1 -1
114,6 20,25
866
260
868
2 -1 1
114,6 69,75
869
197
864
3 1 -1
185,4 20,25 2266
1404
2268
4 1 1
185,4 69,75 2253
2011
2258
5 -1.414 0 100 45
1020
405 1023
6 1.414 0 200 45
4086
4858 4090
7 0
-1.414
150
10
776
121
780
8 0
1.414
150
80
790
423
789
9 0 0
150
45
2316
1406
2301
Na podstawie danych zgromadzonych w tablicy 6.3 obliczamy:
1) wariancję błędów aproksymacji (ze wzoru 6.41):
σ
a
2
= 536,7
2) wariancję wyników pomiarów (ze wzorów 6.42):
σ
2
= 1232
3) wartość funkcji testowej (ze wzoru 6.40): F = 0,4358
Dla rozkładu F – Snedecora odczytujemy z tablicy statystycznej wartość krytyczną dla
poziomu ufności
α = 0,9 i liczby stopni swobody f
1
= n(r - 1) = 81 oraz f
2
= n – N
b
= 6.
Wartość ta wynosi
1
2
,
,
f
f
F
α
= 1,8482.
Wobec spełnienia warunku F ≤
1
2
,
,
f
f
F
α
nie mamy podstaw do odrzucenia hipotezy
zerowej i przyjmujemy, że model obiektu jest adekwatny.
100
120
140
160
180
200
20
40
60
80
0
500
1000
1500
2000
2500
3000
3500
4000
4500
kąt ustawienia działa [deg]
prędkość początkowa [m/s]
zasięg
[m]
Rys. 6.10. Identyfikacja modelu zjawiska ruchu pocisku balistycznego
Optymalizację funkcji modelu przeprowadzono w oparciu o metodę gradientową
i uzyskano następujące wartości zmiennych wejściowych v
0opt
= 200m/s i
ϕ
opt
= 45
°.
22
Maksymalny zasięg pocisku balistycznego obliczony z funkcji modelu wynosi zatem
s
max
= 4090m. Rysunki 6.11a i 6.11b przedstawiają cięcia dwuwymiarowej przestrzeni
zmiennych wejściowych v
0
i
ϕ w ekstremum globalnym.
Przykład 2
Rozważmy eksperyment polegający na identyfikacji modelu procesu chemicznego
osadzania warstw półprzewodnikowych z fazy gazowej podczas produkcji układów
scalonych. Analizujemy wpływ dwóch wielkości wejściowych x
1
i x
2
na jednorodność
warstwy półprzewodnikowej (zmienna wyjściowa y), obliczonej na podstawie pomiarów
rezystancji warstwy za pomocą sondy czteroostrzowej w 49 punktach rozmieszczonych na
powierzchni płytki półprzewodnikowej.
Wielkościami wejściowymi są:
1) ciśnienie w jednostkach [torr] (1 torr oznacza ciśnienie 1 mm słupa rtęci) – zmienna x
1
,
2) iloraz objętości dwóch związków chemicznych wchodzących w skład mieszaniny
gazowej (H
2
i WF
6
) – zmienna x
2
.
Przyjęto następujące granice zmiennych wejściowych:
1) x
1min
= 4 torr, x
1max
= 80 torr,
2) x
2min
= 2, x
2max
= 10.
Celem badań jest wyznaczenie funkcji modelu procesu w postaci wielomianu
algebraicznego drugiego stopnia o postaci:
2
2
22
2
1
11
2
1
12
2
2
1
1
0
ˆ
x
a
x
a
x
x
a
x
a
x
a
a
y
+
+
+
+
+
=
(6.45)
Na podstawie charakterystyki obiektu badań i założonego celu badań wybrano
kompozycyjny, centralny plan doświadczenia o wartości ramienia gwiezdnego
α = 0,71
zawierający n = 11 układów (tablica 6.4). Dla każdego układu wykonano pomiary
rezystancji w 49 punktach i wyznaczono wartości parametrów statystycznych: wartość
średnią
u
y
i odchylenie standardowe
σ
u
. W tablicy 6.4 podano wartości odchylenia
standardowego podzielone przez wartość średnią i wyrażone w procentach.
Rozwiązując układ równań normalnych (6.35) uzyskano następujące wartości
parametrów modelu dla unormowanych zmiennych wejściowych:
a
0
= 6,198 a
1
= -1,910 a
2
= -0,224 a
12
= 1,686 a
11
= -0,197 a
22
= -0,297
10
20
30
40
50
60
70
80
1000
1500
2000
2500
3000
3500
4000
4500
kąt ustawienia działa [deg]
zasięg
[m]
100
110
120
130
140
150
160
170
180
190
200
1000
1500
2000
2500
3000
3500
4000
4500
prędkość początkowa [m/s]
zasięg
[m]
Rys. 6.11. C
ięcia dwuwymiarowej przestrzeni zmiennych wejściowych
w ekstremum globalnym
: a) s = f(
ϕ
) dla
v
0
= v
0opt
, b) s = f(
v
0
) dla
ϕ =ϕ
opt
23
Tablica 6.4.
wielkość
wyjściowa
unormowane
wartości
zmiennych
wejściowych
rzeczywiste
wartości
zmiennych
wejściowych
numer
układu
u
1
u
x(
2
u
x(
x
u1
x
u2
%
100
2
⋅
u
u
y
σ
1 -0,71
-0,71
15,13
3,17 8,6
2 -0,71
0,71
15,13
8,83 6,9
3 0,71
-0,71
68,87
3,17 3,4
4 0,71
0,71
68,87
8,83 5,1
5 -1 0 4 6 7,3
6 1 0 80
6 4,6
7 0 -1
42
2 6,3
8 0 1 42
10 5,4
9 0 0 42
6 6,2
10 0 0 42
6 6,4
11 0 0 42
6 5,0
Następnie przeprowadzono analizę istotności współczynników funkcji aproksymującej
testem t-Studenta. Współczynniki a
11
i a
22
okazały się nieistotne więc dokonano
identyfikacji modelu o uproszczonej postaci:
2
1
12
2
2
1
1
0
ˆ
x
x
a
x
a
x
a
a
y
+
+
+
=
(6.46)
uzyskując następujące wartości parametrów:
a
0
= 5,978 a
1
= -1,910 a
2
= -0,224 a
12
= 1,686
Adekwatność przyjętej struktury modelu zweryfikowano testem F-Snedecora.
Pozytywne przejście tego testu umożliwiło eliminację współczynników a
11
i a
22
z funkcji
modelu.
Ostatecznie korzystając z operacji denormowania uzyskano następujące wartości
współczynników dla rzeczywistych zmiennych wejściowych:
a
0
= 11,246 a
1
= -0,117 a
2
= -0,526 a
12
= 0,011
0
20
40
60
80
2
4
6
8
10
2
3
4
5
6
7
8
9
10
iloraz objętości H
2
/WF
6
ciśnienie [torr]
jednorodność
Rys. 6.12.
Identyfikacja modelu procesu
osadzania warstw półprzewodnikowych