Interpolacja powierzchni odpowiedzi metodą Kriging

W wielu dziedzinach inżynierii najczęściej spotyka się powierzchnie odpowiedzi o charakterze nieliniowym. Dotyczy to zwłaszcza dziedziny montażu elektronicznego. Wynika to z istnienia zjawisk oraz materiałów o właściwościach nieliniowych. Jako model powierzchni odpowiedzi stosuje metody stochastyczne. Przykładem jest metoda Kriging, która należy do metod interpolacyjnych, w związku z czym można ją stosować w przypadku braku czynników szumu.

W przypadku metod stochastycznych przyjmuje się założenie, że uzyskana w wy-niku przeprowadzonego eksperymentu powierzchnia odpowiedzi ma charakter deter-ministyczny y(X), a jednocześnie może być opisana jako realizacja procesu stocha-stycznego. Stochastyczny charakter procesu dotyczy czynnika błędu. Przykładem takiego procesu mogą być obliczenia numeryczne. W takim przypadku model po-wierzchni odpowiedzi można potraktować jako kombinację modelu matematycznego, np. wielomianu oraz dodatkowego czynnika określanego czynnikiem błędu ε(x):



y x=

∑

i=1 k

_i f_ixx . (3.75)

Błąd ε(x) jest odchyleniem wyników eksperymentu od założonego modelu mate-matycznego, co przedstawiono na rysunku 3.32 [219, 220]. Błąd ten można potrak-tować jako wynik realizacji procesu stochastycznego o wartości średniej ε równej 0 i kowariancji V między dwoma punktami przestrzeni u i v:

V u, v=²R u ,v , (3.76) gdzie σ2 – wariancja procesu, R(u,v) – korelacja.

Funkcja kowariancji odpowiada za gładkość dopasowania powierzchni odpowiedzi. W przypadku gładkich powierzchni odpowiedzi funkcja kowariancji ma pochodne, w przeciwieństwie do powierzchni nieregularnych.

Rys. 3.32. Metoda dopasowanie powierzchni odpowiedzi w przypadku metod stochastycznych

Procedura dopasowania modelu powierzchni odpowiedzi składa się z dwóch eta-pów:

● wyznaczenia funkcji głównej np. metodą regresji,

● interpolacji odchyleń w punktach eksperymentu, tak jakby nie została wcześniej przeprowadzona procedura regresji.

Metoda Kriging jest stosunkowo trudna w implementacji. Wymaga ona precyzyjne-go oszacowania parametrów modelu. Z kolei zaletą tej metody jest możliwość oszaco-wania modelu powierzchni błędu odpowiedzi. Powierzchnia błędu pozwala na za-stosowanie procedury sekwencyjnej, tzn. oszacowanie położenia kolejnego punktu w przestrzeni czynnikowej, który zapewnia najlepszą poprawę modelu powierzchni od-powiedzi.

Metoda Kriging zapewnia stosunkowo dużą elastyczność ze względu na dobór od-powiedniej funkcji korelacji R(u,v). Dodatkowo w niektórych przypadkach – w zależ-ności od funkcji korelacji – można uzyskać model powierzchni odpowiedzi na zasadzie aproksymacji. W przypadku eksperymentów numerycznych, które z natury są deter-ministyczne, najczęściej stosuje się jednak metodę interpolacji.

Jeżeli macierz x oznacza punkty wybranego planu eksperymentu dla k czynników oraz n punktów, natomiast macierz y wyniki eksperymentu, to można zapisać, że:

yⁱ=f



x₁ⁱ , x₂ⁱ , , x_kⁱ



; i=1, , n , (3.77) gdzie

x =

∣

x₁1 x₂1 ⋯ x_k1 x₁2 . . . ⋮ . . . x₁ⁿ . . x_kⁿ

∣

, y=

∣

y1 y2 ⋮ yⁿ

∣

. (3.78)

Zgodnie z zależnością 3.75 można dopasować model powierzchni odpowiedzi do wyników doświadczeń według następującej zależności:

yi=

∑

h=1 k

_hf_h



xi



ⁱ, i=1,, n , (3.79) gdzie fh(xi) może być liniową lub nieliniową funkcją położenia xi, βh – nieznany współ-czynnik proporcjonalności założonej funkcji, εi – czynnik błędu odchylenia założonej funkcji od wyników eksperymentu – błąd ten ma rozkład normalny N(0, δ2).

Dodatkowo można przyjąć założenie, że czynniki błędu εi można zapisać jako funk-cję ciągłą położenia xi:

ⁱ=ⁱ



xⁱ



. (3.80)

W takim przypadku czynniki błędu będą wzajemnie skorelowane oraz zależne od wartości odpowiedzi i odległości między punktami w przestrzeni projektowej. Jeżeli punkty te znajdują się blisko siebie, wówczas błędy powinny być podobne, co oznacza dużą wartość korelacji. Z tego też powodu przyjmuje się założenie, że korelacja błę-dów jest funkcją odległości między punktami. Za funkcję odległości można przyjąć zwykłą odległość euklidesową lub częściej następującą funkcję o współczynnikach wa-gowych Θ i p:

d



xi, xj



=

∑

h=1 k

_h

∣

x_hi−x_hj

∣

^ph, (3.81) gdzie współczynnik Θ jest większy od 0, natomiast współczynnika ph zawiera się w przedziale [1, 2].

Korzystając z przyjętej funkcji odległości opisanej zależnością (3.81), korelację między błędami można wyznaczyć z zależności:

corr

[





xⁱ



, 



x^j

]

= 1

e^{dxi, xj . (3.82)} Funkcja korelacji ma następujące właściwości:

● jeżeli odległość między punktami jest mała, to wartość funkcji korelacji jest duża,

● jeżeli odległość między punktami jest duża, to wartość funkcji korelacji zbliża się do 0.

Wartości funkcji korelacji tworzą, tzw. macierz korelacji R o wymiarach n×n, która ma kluczowe znaczenie dla modelu powierzchni odpowiedzi:

R=

∣

r_1,1 ⋯ r_{1, n}

⋮ . .

r_n,1 . r_{n , n}

∣

, (3.83)

gdzie współczynniki ri,j są wektorami postaci:

r_{i , j}=corr

[





xⁱ



, 



x^j

]

. (3.84)

Funkcja korelacji oraz macierzy korelacji R pozwala na wyznaczenie prostego li-niowego modelu powierzchni odpowiedzi, np. metodą regresji lub rachunku wariacyj-nego. Wyznaczony model powierzchni odpowiedzi ma kilka ważnych zalet, z których najważniejsza to możliwość zastąpienia czynników regresji współczynnikami o stałej wartości μ:

=

∑

h=1 k

_hf_h



xⁱ



, i=1,, n . (3.85) Pozwala to, m.in. na przepisanie zależności (3.79) w postaci:

y



xⁱ



=



xⁱ



, 



xⁱ



N



0,²



. (3.86) Metoda Kriging jest stosowana do interpolacji wyników eksperymentu. Niestety jej główną wadą jest implementacja praktyczna, która wynika z następujących proble-mów:

● konieczność estymacji współczynników funkcji odległości opisanej zależnością (3.81); w praktyce dla przypadku jednowymiarowego należy oszacować 2k+2 współczynników, to znaczy: μ, δ2, θ1 ... θk, p1 ... pk,

● poprawna definicji funkcji odległości oraz estymacja jej współczynników w przy-padku wielowymiarowym.

Estymację współczynników funkcji odległości można otrzymać metodą największej wiarygodności przez maksymalizację funkcji F opisanej zależnością:

F = ¹ 2n/ 2²^{n/ 2}

∣

R

∣

^{1/ 2} 1 e y−1 ' R−1y−1 22 , (3.87)

gdzie 1 – n-wymiarowy wektor jednostkowy, y – n-wymiarowy wektor zaobserwowa-nych odpowiedzi:

Estymatory współczynników μ oraz δ2, które maksymalizują funkcję F opisaną za-leżnością (3.87), można zapisać w postaci:

=1' R⁻¹y 1 ' R⁻¹1  ²=^y−1 ' R⁻¹y−1  n . (3.89)

Z kolei najlepszym nieobciążonym estymatorem wartości odpowiedzi y w punkcie x' jest estymator opisany zależnością:

y x' = r ' R⁻¹y−1  , (3.90) gdzie r – n-wymiarowy wektor macierzowym, którego współczynniki można wy-znaczyć z zależności:

r_ix ' =corr

[

x∧, xⁱ

]

. (3.91)

Jak wcześniej podano, zaletą metody Kriging jest możliwość oszacowania błędu dopasowania modelu powierzchni odpowiedzi. Oszacowanie takie można wykonać przez porównanie wartości odpowiedzi y(x') oraz wartości rzeczywistej yi w tym punk-cie, korzystając z definicji wartości średniej błędu kwadratowego s2(x') w postaci:

s2x ' =E

[

y x ' − y x ' 



²

]

= ²

[

1−r ' R r 1^{−1 ' R−1r '}

1' R⁻¹1

]

. (3.92) Dla większości przypadków wygodniej jest jednak posługiwać się pierwiastkiem wartości średniej błędu kwadratowego s(x):

s x=



s x^{2. (3.93)}

Jak wynika z zależności (3.92), oszacowanie błędu interpolacji jest funkcją macierzy korelacji R opisanej zależnościami (3.83) i (3.84). Elementy macierzy R można wyznaczyć dla przypadku jednowymiarowego z zależności:

r_{i , j}= 1

e^h∣xi−xj

h^p∣ . (3.94)

Aby wyznaczyć macierz korelacji R należy dokonać estymacji współczynników θh

oraz ph, co stanowi ważny problem metody Kriging. W tym celu można przyjąć założe-nie, że współczynnik p jest równy 2, natomiast wartość współczynnika Θ wyznacza się wówczas – ze względu na brak odpowiedniego estymatora – wybraną metodą optymalizacji. Współczynnik Θ można uznać za miarę istotności zmiennej xh, co można wyrazić w postaci stwierdzenia, że nawet małe wartości odległości dh^i,j:

d^{i , j}_h =

∣

xⁱ_h−x j_h

∣

, (3.95) mogą prowadzić do dużych zmian wartości funkcji w punktach xi i xj. Wartość współ-czynnika Θh ma istotne znaczenie:

● jeżeli wartość Θh jest duża, to małe zmiany dh^i,j prowadzą do małej wartości korela-cji,

● jeżeli wartość Θh jest mała, to małe zmiany dh^i,j prowadzą do dużej wartości korela-cji,

W praktyce stosuje się następujące metody estymacji współczynnika Θh:

● metoda oparta na procedurze sprawdzianu krzyżowego [123, 221],

● metoda oparta na wyznaczniku macierzy korelacji det R [222],

● metoda zaproponowana przez autora, oparta na rachunku wariacyjnym [208, 206]; metoda polega na wyznaczeniu wartości maksymalnej funkcjonału parametryczne-go F(Θ,p) opisująceparametryczne-go całkowity błąd interpolacji.

Metoda oparta na procedurze sprawdzianu krzyżowego (ang. cross-validation) po-lega na wyznaczeniu błędu dopasowania uzyskanej powierzchni odpowiedzi, według schematu: E  x=



1 n

∑

i=1 n



y_−ix − y x 



², (3.96)

gdzie n – liczba punktów eksperymentu, y x – dopasowana powierzchnią od-powiedzi, y_−ix – dopasowana powierzchnię powierzchni uzyskaną dla n – 1 punk-tów, przy założeniu, że „usunięty” jest i-ty punkt eksperymentu.

Procedura wymaga zatem dopasowania n – 1 powierzchni odpowiedzi y_−ix  . Optymalną wartość współczynnika Θuzyskuje się w wyniku minimalizacji błędu do-pasowania E(x) opisanego zależnością (3.96), tj.:

 =min

 E  x . (3.97)

Metoda oparta na wyznaczniku macierzy korelacji det R polega na wyznaczeniu wartości minimalnej funkcjonału opisanego zależnością:

F  , p=det R

1

n_2. (3.98)

Jeżeli przyjmie się założenie, że wartość parametru p = 2, to estymator parametru Θ można zapisać w postaci:

 =min  F  , p=2=min  det R 1 n_2. (3.99)

Z kolei metoda zaproponowana przez autora polega na wyznaczeniu wartości mak-symalnej funkcjonału parametrycznego F(Θ, p) w postaci całki oznaczonej opisującej całkowity błąd interpolacji:

F



, p



=

∫

x ' x ' '

s



x ,, p



d x , (3.100)

gdzie x – wektor punktów, w których szacuje się wartość odpowiedzi, s(x) – pierwiastek z wartości średniej błędu kwadratowego.

W przypadku przyjęcia założenia, że wartość parametru p = 2, estymator parametru Θ można zapisać jako:

 =min  F , p=2=min 

∫

x ' x ' ' s



x ,, 2



d x . (3.101)

Wartość estymatora parametru Θ wyznacza się metodą numeryczną. Metoda ta po-lega na numerycznym całkowaniu oraz optymalizacji z wykorzystaniem metody sim-pleks [183]. W takim przypadku funkcjonał opisany zależnością (3.100) należy przepi-sać w postaci:

F



, p



=

∑

x ' x ' '

s



x , , p



d x (3.102)

Innym istotnym problemem metody Kriging jest estymacja współczynników Θhi

w przypadku wielowymiarowym. Istnieje wówczas konieczność oszacowania parame-trów Θhi dla każdego wymiaru oddzielnie. Dla przypadku dwuwymiarowego funkcję odległości można zapisać w postaci:

r_{i , j}= 1

e^h1∣x₁i−x₁j

h1^p∣_h2∣x₂i−x₂j

h2^p∣ , (3.103)

co wymaga estymacji dwóch współczynników: Θh1 i Θh2.

Problem estymacji wielowymiarowej występuje stosunkowo często w praktyce in-żynierskiej i stanowi poważne wyzwanie dla istniejących obecnie metod. Wydaje się, że jedną z najbardziej optymalnych metod w tym celu, jest metoda sprawdzianu krzy-żowego oraz metoda zaproponowana przez autora. Przykładem odpowiedzi, która wy-maga takiego podejścia może być powierzchnia odpowiedzi, która w jednym wymiarze ma charakter monotoniczny, w drugim natomiast charakter okresowy, np. model funk-cji opisanej następującą zależnością (rys. 3.33):

Rys. 3.33. Przykład modelu powierzchni odpowiedzi, wymagającej estymacji dwóch parametrów Θ

W dokumencie Numeryczne metody projektowania termomechanicznego w montażu elektronicznym (Stron 144-151)