Aapproximation of the retardation spectrum of viscoelastic plant materials using Lagendre functions

Acta Agrophysica, 2005, 6(3), 795-813

ALGORYTM APROKSYMACJI SPEKTRUM RETARDACJI ROŚLINNYCH MATERIAŁÓW LEPKOSPRĘśYSTYCH FUNKCJAMI LEGENDRE’A

Anna Stankiewicz

Katedra Podstaw Techniki, Akademia Rolnicza ul. Doświadczalna 50A, 20-280 Lublin

e-mail: anna.stankiewicz@ar.lublin.pl

S t r e s z c z e n i e . W pracy rozwaŜa się problem wyznaczania spektrum retardacji liniowych mate-riałów lepkospręŜystych na podstawie dyskretnych zakłóconych pomiarów czasowych przebiegów funkcji pełzania (retardacji), zgromadzonych w standardowym teście pełzania. Zaproponowano metodę aproksymacji ciągłego spektrum częstotliwości retardacji skończoną sumą ortonormalnych funkcji Legendre’a. Problem wyznaczania spektrum retardacji jest źle uwarunkowanym problemem odwrot-nym. Do stabilizacji jego rozwiązania zastosowano technikę regularyzacji Tichonowa ze współczyn-nikim regularyzacji dobranym uogólnioną metodą sprawdzania krzyŜowego (GCV). Pokazano, Ŝe dokładność aproksymacji spektrum retardacji zaleŜy zarówno od zakłóceń w pomiarach funkcji pełzania oraz wartości współczynnika regularyzacji, jak i odpowiedniego doboru parametrów i liczby funkcji bazowych. Załączono wyniki eksperymentu numerycznynego.

S ł o w a k l u c z o w e : lepkospręŜystość, spektrum retardacji, algorytm identyfikacji, regularyzacja, funkcje Legendre’a

WSTĘP

Spektrum retardacji, lub równowaŜnie spektrum relaksacji, fundamentalne dla modeli konstytutywnych ośrodków liniowo lepkospręŜystych, niesie pełną infor-mację o właściwościach mechanicznych tych materiałów. Spektra relaksacji i retar-dacji są uŜytecznym narzędziem badania materiałów lepkospręŜystych, poniewaŜ na ich podstawie moŜna wyznaczyć dowolne inne liniowe charakterystyki opisujące materiały w dziedzinie czasu a takŜe w dziedzinie częstotliwości, w szczególności te charakterystyki, które są wykorzystywane współcześnie w obliczeniach inŜynier-skich, wymagających jak FEM i BEM szczegółowej wiedzy o zjawiskach dyna-micznych zachodzących w materiałach [11,14].

Zarówno spektrum retardacji jak i spektrum relaksacji nie są wprost dostępne pomiarowo, muszą więc być wyznaczane w oparciu o pomiarowo dostępne charakterystyki materiałów. Większość znanych w literaturze metod wyznaczania spektrum relaksacji lub retardacji, opracowanych głównie dla polimerów, bazuje na danych zgromadzonych w wyniku testu, w którym materiał poddawany jest odkształceniom oscylacyjnym o małej amplitudzie w dostatecznie duŜym zakresie częstotliwości [4,12,14-15]. Tylko kilka metod, przeznaczonych głównie do wyznaczania spektrum relaksacji materiałów biologicznych, opartych jest o po-miary przeprowadzane w dziedzinie czasu [5,8,16,18-20,21]. Metody Ter Haara [21] i Fujihary i in. [5] nie gwarantują wyznaczenia spektrum relaksacji z zada-walającą dokładnością, przede wszystkich jednak nie uwzględniają złego uwarun-kowania zadania identyfikacji spektrum, co praktycznie uniemoŜliwia ich zastosowanie w przypadku zakłóconych danych pomiarowych. W pracach [18-20] przedstawiono algorytmy identyfikacji spektrum czasów i częstotliwości relaksacji, w których dla zagwarantowania dobrego uwarunkowania zadania zastosowano róŜne techniki stabilizacji. W literaturze brak jest natomiast algorytmów iden-tyfikacji spektrum retardacji na podstawie danych zgromadzonych w dziedzinie czasu w teście pełzania, który jest dla materiałów pochodzenia roślinnego podsta-wowym źródłem informacji o ich właściwościach mechanicznych [2,17].

W tej pracy zaprezentowano metodę wyznaczania spektrum częstotliwości retardacji na podstawie dyskretnych pomiarów funkcji pełzania zgromadzonych w standardowym teście pełzania, bazującą na sprowadzeniu ciągłego zagadnienia identyfikacji spektrum retardacji do liniowo-kwadratowego problemu dyskretnego za pomocą rozwinięcia spektrum względem ortonormalnej bazy w przestrzeni funkcyjnej. Jest to podejście stosowane zarówno w teorii aproksymacji, jak i w za-daniach identyfikacji obiektów dynamicznych z czasem ciągłym [1].

MATERIAŁ I METODY Spektrum retardacji

W zakresie niewielkich deformacji związek między odkształceniem ε

( )

t a

naprę-Ŝeniem σ

( )

t w materiale lepkospręŜystym opisuje całkowe równanie konstytutywne oparte o zasadę superpozycji Boltzmanna [3,14]

( )

(

τ

) ( )

σ τ τ ε t J t d t &

∫

∞ − − =

gdzie J

( )

t jest funkcją pełzania (retardacji). Dla lepkospręŜystych materiałów po-chodzenia roślinnego, takich jak wysoko uwodnione warzywa i owoce, funkcja pełzania J

( )

t jest dana modelem Kelvina [2,17]

( )

t J L

( )

ν

(

e ν

)

dν

J = +∞

∫

− −t

0

0 1 (1)

gdzie J₀ to podatność natychmiastowa, natomiast funkcja L

( )

ν oznacza spektrum częstotliwości retardacji, które charakteryzuje frakcję elementów modelu (1) o często-tliwościach retardacji zawierających się pomiędzy ν a ν +dν [3,14].

Identyfikacja [1] spektrum retardacji L

( )

ν sprowadza się do numerycznego problemu wyznaczenia rozwiązania całkowego równania Fredholma 1-go rodzaju (1) na podstawie dyskretnych danych pomiarowych. Praktyczna trudność w iden-tyfikacji spektrum retardacji ma źródło w problemie matematycznym, bowiem jak wiadomo zadanie to jest źle postawionym problemem odwrotnym [14,15,22] i jego rozwiązanie wymaga stosowania specjalnych metod. W tej pracy zapropo-nowano metodę identyfikacji spektrum retardacji, w której spektrum retardacji

( )

ν

L przybliŜa się kombinacją liniową ortonormalnych funkcji Legendre'a.

Podstawy matematyczne algorytmu

Będziemy zakładać, Ŝe funkcja L

( )

ν ∈L₂

(

0,∞

)

. Funkcje Legendre’a dane wzorem [10]

( )

(

k

)

e P

(

e

)

, k ,, ,

p_k ν = ₂α ₂ +₁ −αν _k1−₂ −2αν =₀₁K ₍₂₎

gdzie P_k

( )

x są wielomianami Legendre’a zdefiniowanymi róŜniczkowym wzorem Rodriguesa [10]

( )

( )

x , k , , , dx d k x P k k k k k 1 01K 2 1 2_{− =} = ! ,

tworzą zupełny ortonormalny układ bazowy w przestrzeni L 0₂

[

,∞

)

funkcji całkowalnych z kwadratem w przedziale

[

0,∞

)

.

Spektrum retardacji L

( )

ν będziemy przybliŜać za pomocą kombinacji liniowej K pierwszych elementów ciągu funkcyjnego

{

p_k

( )

ν

}

postaci

( )

∑

−

( )

= = 1 0 K k k k K g p L ν ν (3)

Oznacza to przybliŜenie funkcji pełzania J

( )

t funkcją

( )

∑

−

( )

= + = 1 0 0 K k k k K t J g t J ϕ (4)

gdzie ϕ_k

( )

t p_k

( )

ν

(

e−tν

)

dν ∞ − =

∫

1 0 (5).

Podstawowe znaczenie dla konstrukcji algorytmu identyfikacji spektrum retardacji ma następujące twierdzenie, które precyzuje postać funkcji bazowych

( )

t

k

ϕ (5). Dowód twierdzenia podano w Dodatku.

Twierdzenie. Niech k ≥0, α >0 i t≥0 . Dla funkcji bazowych Legendre’a

( )

ν k

p funkcje ϕ_k

( )

t (5) dane są wzorem

( ) ( ) ( ) ( ) [ ] ( ) [ i t] , k ,, , t i k k t k i k i k 01K 1 2 1 2 1 2 1 1 2 2 0 1 0 ₌             + + − + + − + = ∏ ∏ = − = α α α α ϕ (6)

lub, równowaŜną, prostą formułą rekurencyjną

( )

( )

_[

[

₍

(

)

₎

]

_]

(

)

(

k

) (

[

(

k

)

)

t

]

, k ,, , k t k t k t k k k t t _k k K 1 0 3 2 1 2 1 3 2 2 4 3 2 1 2 1 2 3 2 1 = + + + + + + + + − + + + = + α α α α ϕ ϕ (7) zaczynając od

( )

(

t

)

t t + = α α ϕ0 2 (8)

Rys. 1. Funkcje bazowe ϕk( )t (6) dla k=0÷3, parametr α = 0,03 i α = 0,3 Fig. 1. The basis functions ϕ_k( )t (6) for k=0÷3, parameter α = 0.03 and α = 0.3

2 0 5 10 1 Czas t Time t, s F u n k cj e b az o w e ϕk (t ), P a -1 B as is f u n ct io n s ϕk (t ), P a -1 α = 0,3 ( )t 0 ϕ ( )t 1 ϕ ϕϕ23( )( )tt 0 0.5 1 2 Czas t Time t, s F u n k cj e b az o w e ϕk (t ), P a -1 B as is f u n ct io n s ϕk (t ), P a -1 α _{= 0,03} ⋅⋅⋅⋅0.5 ϕ0( )t ( )t 1 ϕ ϕϕ23( )( )tt

Funkcje

ϕ

k

( )

t

(7) to funkcje wymierne zmiennej rzeczywistej t ≥0. Przebieg

kilku pierwszych funkcji

ϕ

k

( )

t

dla dwu róŜnych parametrów α przedstawiono

na rysunku 1.

Problem przybliŜania spektrum retardacji L

( )

ν kombinacją liniową funkcji Legendre’a L_K

( )

ν został więc sprowadzony do problemu aproksymacji funkcji pełzania J

( )

t skończoną sumą J_K

( )

t (4) funkcji wymiernych ϕ_k

( )

t (6).

Problem identyfikacji

ZałóŜmy, Ŝe przeprowadzono skończony eksperyment dyskretny, którego rezultatem jest zbiór pomiarów funkcji pełzania J

( )

t_i = J

( ) ( )

t_i +z t_i w chwilach czasu t_i ≥0 , i=1,K,N, gdzie z

( )

t_i jest addytywnym błędem pomiarowym. Jako miarę dokładności modelu (4) będziemy przyjmować klasyczny kwadratowy wskaźnik jakości, który definiując macierz i wektor

( ) ( ) ( ) ( ) ( ) ( )          = − − 1 1 1 1 0 1 1 1 1 1 0 N K N N K N,K t t t t t t ϕ ϕ ϕ ϕ ϕ ϕ K M O M M M K Ψ Ψ Ψ Ψ ( ) ( )         = N N t J t J M 1 J (9)

moŜna zapisać w postaci

( )

[

( )

( )

]

2 1 2 K K N N N i i K i K N J t J t Q g = ∑ − = J −ΨΨΨΨ _, g = (10)

gdzie _gK =

[

g₀ K g_K−1 J₀

]

T _jest

(

_K+₁

)

_{elementowym wektorem}

niezna-nych współczynników w modelu (4), zaś

⋅

oznacza normę euklidesową w przestrzeni RN (dalej takŜe w przestrzeniach RK oraz RK+1).

Problem identyfikacji spektrum retardacji w klasie funkcji L_K

( )

ν postaci (3) sprowadza się więc do rozwiązania, względem g_K , klasycznego liniowego problemu najmniejszej sumy kwadratów

2 1 N N K K RK K min J _, g g Ψ ΨΨ Ψ − + ∈ (11)

Problem identyfikacji ciągłego spektrum retardacji został więc sprowadzony do dyskretnego problemu najmniejszych kwadratów (11). Jak wiadomo [9,18-20, 22] zadanie to jest źle postawione w sensie Hadamarda [22] i nawet bardzo małe błędy w pomiarach mogą spowodować znaczne fluktuacje rozwiązania. Tradycyj-nie stosowanym remedium jest ograniczeTradycyj-nie zbioru rozwiązań dopuszczalnych lub odpowiednia regularyzacja problemu [22]. W tej pracy stosujemy obie

techniki równocześnie. Funkcję L

( )

ν przybliŜono kombinacją liniową funkcji Legendre’a (3), dla zapewnienia dobrego postawienia zadania zastosujemy technikę regularyzacji Tichonowa [22], polegającą na stabilizacji rozwiązania zadania (11) poprzez minimalizację zmodyfikowanego wskaźnika jakości

K K K N N RK K min J _, g g g 2 2 1 − +λ + ∈ ΨΨΨΨ (12)

gdzie λ >0 jest parametrem regularyzacji. Zadanie (12) jest dobrze uwarunko-wane i posiada jednoznaczne rozwiązanie:

(

N KT N K K ,K

)

N KT N

K I J

gλ =ΨΨΨΨ _, ΨΨΨΨ _, +λ _{+1 +1}−1ΨΨΨΨ _, (13) gdzie I_K,K oznacza K×K wymiarową macierzą jednostkową. Wektor gλ_K (13) jest ciągłą funkcją zarówno macierzy ΨΨΨΨ_{N ,K} jak i wektora J_N.

Wzór (13) ma znaczenie głównie teoretyczne, do konstrukcji numerycznego algorytmu wyznaczania macierzy odwrotnej w (13) zastosujemy technikę dekom-pozycji macierzy względem wartości szczególnych (SVD) [9]. Zastosowanie tej techniki uprości takŜe rozwiązanie zadania doboru współczynnika regularyzacji.

Podstawy algebraiczne algorytmu

Niech rozkład SVD macierzy ΨΨΨΨ_N,K przyjmuje postać [9] T K N, UΣΣΣΣ V Ψ Ψ Ψ Ψ = (14)

gdzie macierz ΣΣΣΣ ∈RN,K+1 jest macierzą utworzoną z niezerowych wartości szczególnych σ₁,K,σ_r macierzy ΨΨΨΨ_N,K [9]

(

₁,K, ,0,K,0

)

r diag σ σ = ΣΣΣΣ

(

N,K

)

r =rządΨΨΨΨ , zaś V ∈RK+1,K+1 i U∈RN,N są macierzami ortonormalnymi. Na podstawie wzoru (14), wykorzystując ortonormalność macierzy V , mamy

(

)

(

)

(

)

T ,K K T ,K K T T ,K K K N T K N V I V I V V I , , 1 1 1 1 1 1 + + + + + + + = + = + λ λ λ ΣΣΣΣ ΣΣΣΣ ΣΣΣΣ ΣΣΣΣ Ψ Ψ Ψ Ψ Ψ Ψ Ψ Ψ (15)

a stąd ,wykorzystując ponownie ortonormalność macierzy V , otrzymujemy

gdzie ΛΛΛΛ jest

(

K+1

) (

× K +1

)

wymiarową macierzą diagonalną postaci

(

)

(

)

(

1 σ₁2+λ ,K,1 σ2+λ ,1λ,K,1λ

)

=diag r Λ ΛΛ Λ (17)

Ostatecznie więc na podstawie wzorów (13) oraz (15)-(17) parametr gλ_K dany jest formułą

gλ_K =VΛΛΛΛVTΨΨΨΨ_N,KT J_N (18)

Dobór parametru regularyzacji

Skuteczność stosowania techniki regularyzacji Tichonowa zaleŜy od metody doboru parametru regularyzacji. Dobór współczynnika regularyzacji jest sztuką wymagającą zarówno intuicji, dobrych heurystyk, jak i często wiedzy aprio-rycznej o zakłóceniach. Generalnie parametr λ dobiera się tak, aby uzyskane rozwiązanie gλ_K było moŜliwie bliskie rozwiązaniu, jakie uzyskalibyśmy dla pomiarów bezszumowych. Wśród wielu metod zaproponowanych dla doboru parametru regularyzacji λ najczęściej wymienia się: zasadę nie zrównowaŜenia Morozowa, technikę L-krzywej, kryterium quasi-optymalności oraz uogólnioną metodę sprawdzania krzyŜowego (Generalized Cross Validation-GCV). Dwie pierwsze metody wymagają znajomości wariancji zakłóceń obciąŜających pomiary, w pozostałych nie jest to konieczne. W pracy [18] zastosowano prostą regułę doboru współczynnika regularyzacji, która uwzględniając wprost wartość wskaź-nika jakości identyfikacji, gwarantuje załoŜoną dokładność modelu, podczas gdy w [19] parametr λ dobiera się tak, aby zagwarantować z góry załoŜoną wartość normy z gλ_K.

W tej pracy zastosujemy technikę GCV, o jej wyborze decyduje zarówno prosta idea, elegancja matematyczna oraz łatwy w implementacji algorytm jak I to, Ŝe dla zastosowania GCV nie jest niezbędna Ŝadna wiedza aprioryczna o zakłóceniach w pomiarze funkcji pełzania. Schemat doboru współczynnika regularyzacji metodą sprawdzania krzyŜowego (CV) został zaproponowany w [23], a następnie uogólniony w [7] i [13]. Idea CV jest następująca: dla danego parametru regularyzacji λ dzielimy zbiór danych na dwa podzbiory, jeden z nich jest wykorzystywany do wyznaczenia rozwiązania zregularyzowanego, drugi podzbiór zaś wykorzystywany jest do oceny błędu aproksymacji. Parametr λ dobiera się tak, aby minimalizował on sumę błędów dla danych z drugiego podzbioru dla wszystkich moŜliwych podziałów. W uogólnionej zasadzie spraw-dzania krzyŜowego, zachowując ideę CV, dane poddaje się pewnemu przekształ-ceniu unitarnemu bazującemu na macierzy Fouriera [7,13], co oznacza, Ŝe GCV jest wersją CV inwariantną ze względu na obrót danych o wielokrotności kąta

N

π

2 . W pracach [7,13] pokazano, Ŝe dla układów liniowych parametr regulary-zacji zgodnie z zasadą GCV dobiera się jako globalne minimum

( )

λ λ λ GCV GCV arg minV 0 > = (19)

tzw. funkcji GCV danej wzorem

( )

_λ

( )

_λ 2

(

_trace

[

_M

( )

_λ

]

)

2

V_GCV = r gdzie M

( )

λ jest macierzą postaci

( )

(

)

T K N ,K K K N T K N K N N N, M λ = I −ΨΨΨΨ _, ΨΨΨΨ _, ΨΨΨΨ _, +λ I _{+1 +1}−1ΨΨΨΨ _,

natomiast r

( )

λ = M

( )

λ J_N = J_N−ΨΨΨΨ_N,K g_Kλ jest wektorem residualnym dla rozwiązania zregularyzowanego (13).

Wykorzystując rozkład SVD macierzy ΨΨΨΨ_N,K (14), oraz wzory (16)-(18) i definiując N wymiarowy wektor Y =UT J_N, występującą w liczniku funkcji

( )

λ GCV

V normę z wektora residualnego r

( )

λ moŜna łatwo obliczyć jako funkcję wartości szczególnych σ₁,K,σ_r oraz składowych wektora

Y

postaci:

( )

(

)

∑

∑

+ = = ₊ + = N r i i r i i i _y y 1 2 1 2 2 2 2 2 λ σ λ λ r

RównieŜ ślad występujący w mianowniku funkcji V_GCV

( )

λ dany jest prostym wzorem

( )

[

]

∑ ∑ = = + = − + + − = r i _i r i _i i _{N r} N M trace 1 2 1 2 2 λ σ λ λ σ σ λ

a stąd funkcja GCV dana jest dogodnym wzorem analitycznym

( )

(

)

2 1 2 1 2 1 2 2 2 2         + + −         + + =

∑

∑

∑

= + = = r i _i N r i i r i i i GCV y N r y V λ σ λ λ σ λ λ (20)

W pracy [7] pokazano, dla przypadku probabilistycznego, Ŝe funkcja

( )

λ GCV

EV posiada w przedziale

[

0,∞

)

przynajmniej jedno minimum lokalne 0

≥

∗

λ . Generalnie jednak problem istnienia rozwiązania zadania optymalizacji (19) nie został rozstrzygnięty ani w pracach twórców metody i ich współpra-cowników [8,13], ani w licznych pracach poświęconych zastosowaniom zasady

minimum właściwe w przedziale

[

0,∞

)

jest istnienie 1≤i≤r, takiego Ŝe 0

≠

i

y , w przeciwnym bowiem przypadku V_GCV

( )

λ (20) jest funkcją monoto-nicznie malejącą dla λ∈

[

0,∞

)

. PoniŜej wykazano, iŜ dla rozpatrywanego zadania identyfikacji spektrum retardacji warunek ten jest spełniony.

Własność 1. Niech K ≥1, N ≥1 i α >0 . Wówczas istnieje 1≤i≤r, takie Ŝe 0

≠

i

y .

D o w ó d przeprowadzimy nie wprost. ZałóŜmy, Ŝe y_i =0 dla i=₁,K,r_{. Niech}

N T T N T _J U U J U Y Y Y       = = =       2 1 2 1 (21)

gdzie Y₁∈Rr, Y₂∈RN−r, U₁∈RN,r, a U₂∈RN,N−r będzie rozkładem blokowym iloczynu Y =UT J_N. RównieŜ

rozkład SVD

(14)

macierzy

ΨΨΨΨ_{N ,K}

moŜna przedstawić w postaci blokowej

[

]

            = = T T T K N 2 1 11 2 1 V V U U V U , 0 0 0 ΣΣΣΣ ΣΣΣΣ Ψ ΨΨ Ψ (22)

gdzie macierz ΣΣΣΣ₁₁=diag

(

σ₁,K,σ_r

)

, zaś V₁∈RK+1,r a V₂∈RK+1,K+1−r. Wyznaczymy iloczyn Y V J U V J W =ΨΨΨΨ_N,KT _N = ΣΣΣΣT T _N = ΣΣΣΣT

Na podstawie (22) i (21) wykorzystując załoŜenie Y₁=0_r, gdzie 0_r oznacza

r wymiarowy wektor zerowy, otrzymujemy

1 1 11 1 1 11 1 = = + = = T T _K K N J V U J V Y 0 W ΨΨΨΨ _{, N} ΣΣΣΣ _N ΣΣΣΣ (23)

Równocześnie uwzględniając definicje macierzy ΨΨΨΨ_{N ,K} (9) i wektora J_N (9), łatwo sprawdzić, Ŝe W_K+1 =∑_iN₌₁J

( )

t_i , a poniewaŜ J

( )

t_i >0 dla i=₁,K,N_,

równieŜ W_K+1 >0, co jest sprzeczne z (23) i kończy dowód.

□

Łatwo zauwaŜyć, iŜ dowód tej własności wykorzystuje zarówno strukturę macierzy ΨΨΨΨ_{N ,K} o jednostkowych elementach ostatniej kolumny, która wynika ze struktury modelu Kelvina (1), jak i specyfikę procesu pełzania, w którym

( )

t_i > 0

J . Prosty wystarczający warunek a posteriori istnienia rozwiązania zadania (19) podaje kolejna własność.

Własność 2. Niech K ≥1, N ≥1, α >0 i r =rząd

(

ΨΨΨΨ_N,K

)

< N. Jeśli

[

]

[ ]

2 1 2 2 1 2 _{N r y N} y N r i i N i i      ∑ ≥ −       ∑ + = = (24)

to istnieje rozwiązanie zadania (19) doboru współczynnika regularyzacji metodą

GCV.

D o w ó d : Funkcja V_GCV

( )

λ jest róŜniczkowalna dla dowolnego λ ≥0, na podstawie wzoru (20) łatwo wyznaczyć formułę analityczną dla pochodnej

( )

λ dλ

dV_GCV , i sprawdzić, Ŝe w szczególności

( )

_[

_]

3 1 2 1 2 ₁ 2 0 r N y d dV r i i N r i i GCV ₋       ∑       ∑ − = = + = σ λ

czyli

dV

_GCV

( )

0

d

λ

<

0

. Jeśli spełniona jest nierówność (24), to

( )

[ ]

[

]

2

( )

0 1 2 2 1 2 GCV N r i i N i i GCV y N y N r V V lim _ − =      ∑ ≥       ∑ = + = = ∞ → λ λ

a to wraz z ciągłością funkcji V_GCV

( )

λ implikuje istnienie jej minimum właści-wego λ_GCV w przedziale

[

0,∞

)

.

Funkcja V_GCV

( )

λ jest róŜniczkowalna dla kaŜdego λ, do rozwiązania zada-nia minimalizacji (19) moŜna więc zastosować dowolną gradientową technikę optymalizacji.

WYNIKI I DYSKUSJA Analiza

Zastąpienie wyjściowego zadania najmniejszych kwadratów (11) zmodyfiko-wanym problemem optymalizacji (12) słuŜyło ograniczeniu fluktuacji rozwiąza-nia gλ_K. Ocenę skuteczności tego podejścia umoŜliwia następujące oszacowanie wynikające wprost z Własności 1 w [18]

(

)

1 2 2 2 1 2 2 2 2 2 _N K r i i i r i GCV i i i K y y GCV _g g < = + =

∑

∑

= = _{σ λ} σ σ λ ₍₂₅₎

gdzie g_KN jest rozwiązaniem normalnym zadania (11). Równość we wzorze (25) ilustruje mechanizm stabilizacji rozwiązania GCV

K

λ

przyj-muje parametr λ_GCV, tym silniej ograniczone są fluktuacje wektora GCV K

λ

g .

Szczegółowe omówienie procesu stabilizacji rozwiązania GCV K

λ

g podano w [18]. PoniewaŜ funkcje bazowe

{

p_k

( )

ν

}

są ortonormalne w przestrzeni L 0₂

[

,∞

)

, dla dowolnego L_K

( )

ν postaci (3) zachodzi równość

( ) ( )

= ∑

[ ]

∑ ∑ = − = − = − = ∞

∫

1 0 2 1 0 1 0 0 2 1 K k k K k K j k j k j K g g p p d g L ν ν ν

gdzie ⋅ ₁ oznacza normę kwadratową w przestrzeni L 0₂

[

,∞

)

. Stąd, poniewaŜ

[

]

T

K K = g0 K g −1 J0

g , otrzymujemy proste oszacowanie

[ ]

[ ]

2 2 0 2 1 0 2 2 1 K K K k k K g J L =

∑

− = g − ≤ g = (26)

Wygładzenie rozwiązania problemu dyskretnego (11) gwarantuje więc równo-cześnie ograniczenie fluktuacji spektrum retardacji L_K

( )

ν , w szczególności zape-wnia ograniczenie fluktuacji wyznaczonego spektrum

( )

∑

−

( )

= = 1 0 K k k K K g p Lˆ ν λGCV ν ₍₂₇₎

Funkcja Lˆ_K

( )

ν (27) jest więc przybliŜeniem w klasie funkcji (3) rzeczywistego spektrum retardacji L

( )

ν optymalnym w sensie kwadratowego wskaźnika jakości identyfikacji Q_N

( )

g_K (10) o ograniczonej fluktuacji Lˆ_K < g_KN

1 .

Spektrum Lˆ_K

( )

ν (27) stanowi oczywiście jedynie przybliŜenie spektrum, które moŜna by uzyskać w klasie funkcji postaci (3), minimalizując wprost (bez regularyzacji) kwadratowy wskaźnik jakości identyfikacji dla pomiarów dokład-nych, czyli przybliŜenie funkcji

( )

∑

−

( )

= = 1 0 K k k N k N K g p L ν ν (28)

gdzie g_KN jest rozwiązaniem normalnym wyjściowego zadania najmniejszych kwadratów (11), jakie uzyskalibyśmy dla pomiarów bezszumowych. PokaŜemy,

Ŝe dokładność aproksymacji zaleŜy zarówno od zakłóceń w pomiarach, jak I współczynnika regularyzacji λ_GCV.

Własność 3. Niech K ≥1, N ≥1 i α >0 . Zachodzą oszacowania N r i _i i GCV N K K N K K y L Lˆ _g GCV _g γ _z σ λ λ _{− ≤ ∑ +} ≤ − =1 3 1 (29) gdzie _zN =

[

z

( )

t₁ K z

( )

t_N

]

T_{, a 1 0} 1 > ∑ = r₌ i σi γ .

D o w ó d : PoniewaŜ g_KN =ΨΨΨΨ_N,K+ J_N, gdzie ΨΨΨΨ_{N ,K}+oznacza pseudoodwrotność Moore’a-Penrose’a [9] macierzy ΨΨΨΨ_{N ,K}, łatwo sprawdzić, Ŝe

(

)

2 1 3 2 1 6 2 2 1 2 2 2 2 2 2         ∑ ≤ ∑ ≤ ∑ + = − = = = r i _i i GCV r i _i i GCV r i GCV i i i GCV N K K y y y GCV σ λ σ λ λ σ σ λ λ _g g (30)

a stąd na podstawie nierówności trójkąta, oszacowania (26) oraz nierówności

N N K N K N K N K g g g z g − ≤ − ≤γ

otrzymujemy następujący ciąg oszacowań

r _N i _i i GCV N K N K N K K N K K N K K y L Lˆ _g GCV _{g g} GCV _{g g g} γ _z σ λ λ λ _{− ≤ − + − ≤ ∑ +} = − =1 3 1 czyli nierówność (29).

□

Na podstawie nierówności (29) wektor współczynników GCV Kλ

g dąŜy do wektora normalnego g_KN, liniowo ze względu na z_N , gdy równocześnie

0

→

GCV

λ i z_N →0. Nierówność (29) gwarantuje równieŜ zbieŜność spektrum retardacji Lˆ_K

( )

ν do LN_K

( )

ν w kaŜdym punkcie ν ich ciągłości, gdy

0

→

GCV

λ i z_N →0.

Z nierówności (29), równości we wzorze (30) oraz definicji stałej γ wynika,

Ŝe dokładność aproksymacji rzeczywistego spektrum retardacji zaleŜy od zakłóceń w pomiarach funkcji pełzania, wartości parametru λ_GCV oraz wartości szczególnych σ₁,K,σ_r macierzy ΨΨΨΨ_K,K, które z kolei zaleŜą zarówno od doboru liczby funkcji bazowych p₀

( )

ν ,K,p_K−1

( )

ν _{, jak i parametru} α_.

Algorytm identyfikacji spektrum retardacji

1. Przeprowadź eksperyment (test pełzania [2,17]) i zgromadź pomiary J

( )

t_i

funkcji pełzania w chwilach czasu t_i ≥0 , i=1, K,N.

3. Wyznacz funkcję V_GCV

( )

λ (20), a następnie wyznacz parametr λ_GCV rozwiązując zadanie optymalizacji (19).

4. Wyznacz macierz

(

K +1

) (

× K+1

)

wymiarową macierz diagonalną

(

)

(

)

(

GCV r GCV GCV GCV

)

diag1 σ₁2+λ ,K,1 σ2+λ ,1λ ,K,1λ = Λ ΛΛ Λ

5. Wyznacz rozwiązanie zregularyzowane GCV K λ g zgodnie ze wzorem N T K N T K V V J gλ = ΛΛΛΛ ΨΨΨΨ _,

6. Wyznacz spektrum częstotliwości retardacji Lˆ_K

( )

ν dane sumą funkcji Legendre’a (27).

U w a g a 1 : Rozkład SVD macierzy ΨΨΨΨ_K,K jest z numerycznego punktu widzenia głównym etapem algorytmu. Procedury SVD są dostępne w prawie kaŜdym pakiecie obliczeniowym (np. svds(A) i svd(A) w Mathcadzie i svd(A) Matlabie). U w a g a 2 : W przedstawionym algorytmie parametr α >0 umoŜliwia zmianę skali czasu. Zachodzi następująca prawidłowość: im mniejszą wartość przyjmuje parametr α tym krótsze są czasy, czyli większe częstotliwości retardacji. Prawidłowość tę ilustruje rysunek 1.

U w a g a 3 : Dobierając parametr α optymalnie moŜemy zagwarantować lepsze dopasowanie modelu do danych eksperymentalnych. W tym celu naleŜy procedurę identyfikacji rozszerzyć o nadrzędny moduł odpowiedzialny za dobór najlepszego parametru α. W praktyce wystarczy jednak zgrubna strategia doboru współczyn-nika α bazująca na porównaniu dla róŜnych wartości stałej α uzyskanych ekspery-mentalnie przebiegów czasowych funkcji pełzania J

( )

t oraz kilku pierwszych elementów ciągu funkcyjnego

{

ϕ_k

( )

t

}

. W ten sposób moŜna takŜe wstępnie oszacować liczbę K składników sumy (3), lub równowaŜnie (4).

Przykład numeryczny

RozwaŜmy materiał lepkospręŜysty, którego spektrum retardacji dane jest rozkładem Gaussa

( )

( )2

( )

2 o 2 2 1 _{ν ν ρ} ρ π ν = e− − L

z parametrami νo =1 i ρ =0.2, podatność natychmiastowa J0 =0.2Pa−1. Przebieg funkcji pełzania J

( )

t = J

( ) ( )

t +z t z addytywnymi zakłóceniami pomiarowymi z

( )

t o rozkładzie jednostajnym w przedziale

[

−0.02,0.02

]

przedstawia rysunek 2(a). Funkcję J

( )

t spróbkowano w N =100 punktach pomiarowych ze stałym okresem próbkowania ∆t=0.06. Parametry K =6 i

α

=0.3 zostały dobrane zgodnie z opisaną powyŜej procedurą, współczynnik regularyzacji λGCV = 0.112⋅10

-2

. Przebiegi „rzeczywistego” spektrum retardacji

( )

ν

L oraz jego przybliŜenia Lˆ_K

( )

ν przedstawiono na rysunku 2(b). Normy z rozwiązania normalnego i zregularyzowanego oraz odpowiednie wartości wskaźnika jakości identyfikacji są następujące:

34 16 gKN = . , gKGCV =1.352 λ _,

( )

2 10 172 1 g = . ⋅ − Q_{N K}N i

( )

g GCV =1.181⋅10−2 K N Q λ .

Rys. 2. Przebieg funkcji pełzania J( )t_i =J( ) ( )t_i +z t_i z zakłóceniami pomiarowymi (a), Spektrum retardacji L( )ν (linia przerywana) i przybliŜenie Lˆ_K( )ν (linia ciągła) (b)

Fig. 2. Creep compliance J( )t_i =J( ) ( )t_i +zt_i corrupted by additive noise (a), Retardation

spectrum L( )ν (dash line) and approximation Lˆ_K( )ν (solid line) (b)

WNIOSKI

1. PrzybliŜenie funkcji L

( )

ν wielomianem funkcyjnym L_K

( )

ν ogranicza zbiór rozwiązań dopuszczalnych do skończenie wymiarowej podprzestrzeni

{

p₀,p₁, , p ₋₁

}

⊂ L₂

(

0,∞

)

Lin K _{K . W konsekwencji problem identyfikacji ciągłego}

spektrum retardacji został sprowadzony do statycznego dyskretnego źle uwarunko-wanego problemu najmniejszych kwadratów. Stabilizację jego rozwiązania zapew-nia zastosowanie techniki regularyzacji Tichonowa, dobór współczynnika regulary-zacji umoŜliwia efektywna obliczeniowo uogólniona metoda sprawdzania krzyŜowego. 2. Dobór ortonormalnych funkcji bazowych gwarantuje, iŜ wygładzenie zregularyzowanego rozwiązania tego problemu zapewnia wygładzenie wyznaczo-nego spektrum retardacji.

Czas t Time t, s 0 5 1 ( )ti J , Pa-1 (a) (b) ( )ν K Lˆ ( )ν L Częstotliwośćν Frequencyν, s-1 0 0 0 1 2 3 0 1 2

3. Wybór parametru α oraz liczby K funkcji bazowych zaleŜy zarówno od wiedzy apriorycznej o badanym procesie, jak i od informacji, którą dysponujemy a

posteriori po zgromadzeniu i wstępnej analizie wyników eksperymentu. Ostate-cznie jednak o wyborze konkretnych wartości tych parametrów decyduje przede wszystkim jakość aproksymacji rzeczywistego spektrum retardacji mierzona wartością wskaźnika jakości identyfikacji, jaką te parametry gwarantują.

4. Zastosowanie funkcji bazowych, dla których całki (5) dane są prostymi formułami analitycznymi, pozwoliło takŜe, co szczególnie istotne w kontekście

źle postawionego problemu odwrotnego, uniknąć błędów przybliŜeń kwadratur numerycznego całkowania, występujących w znanych algorytmach [4,12,15] wyznaczania ciągłego spektrum relaksacji i retardacji.

DODATEK

Dowód Twierdzenia. Łatwo zauwaŜyć, Ŝe całkę ϕ_k

( )

t (5) moŜna przed-stawić jako róŜnicę dwu całek

( )

t p_k

( )

d p_k

( )

e t d _k

( )

_k

( )

t k ν ν ν ν φ φ ϕ =∞

∫

−∞

∫

−ν = 0 − 0 0 (D1)

gdzie funkcja φ_k

( )

t jest zdefiniowana następująco

( )

( )

ν ν φ_k t p_k e−tν d ∞

∫

= 0 (D2) Najpierw pokaŜemy Ŝe całki φ_k

( )

t dane są następującym wzorem

( )

(

)

[

(

)

]

(

)

[

]

K , , k , t i t i k t k i k i k 01 1 2 1 2 1 2 2 0 1 0 ₌ ∏ + + ∏ + − + = = − = α α α φ (D3)

Dowód wzoru (D3) przeprowadzimy metodą indukcji matematycznej. Aby wykazać słuszność wzoru (D3) dla k =0 i k =1 wystarczy, wykorzystując równości P₀

( )

x =1 i P₁

( )

x =x [10], wyznaczyć całki φ₀

( )

t i φ₁

( )

t na podstawie wzoru definicyjnego (D2) i sprawdzić, Ŝe są one dane wzorem (D3). Kolejno

( )

₍

₎

= ∏

[

(

+

)

+

]

+ = ∫ = = − − ∞ 0 0 0 0 2 2 1 1 2 2 i t _{i t} t d e e t α α α α ν α φ αν ν

( )

t = α∞∫

(

− e−

)

e− e−t dν = α

(

α−t

) (

[

α +t

)(

α+t

)

]

φ 6 1 2 αν αν ν 6 3 0 2 1 czyli

( )

(

)

∏

[

(

)

]

∏

[

(

)

]

= = + − + + + = 1 0 0 0 1 2 2 1 2 1 2 1 i i t i t i t α α α φ

ZałóŜmy teraz, Ŝe równość (D3) jest prawdziwa dla

(

k −1

)

i k , gdzie k ≥1. PokaŜemy, Ŝe wzór (D3) jest prawdziwy równieŜ dla

(

k+1

)

Wykorzystując toŜsamość

(

k +1

)

P_k+1

( ) (

x = 2k +1

)

xP_k

( )

x −kP_k−1

( )

x [10, wzór (23,1)] oraz wzory definicyjne (2) i (D2) otrzymujemy równość

(

)

( )

( )

(

)

( )

t k k k t k k t k k t k k k k k 1 1 1 2 3 2 2 1 2 3 2 2 1 2 3 2 1 − + − + − + + + − + + = + φ α φ φ φ czyli równowaŜnie

( )

₍

₎

[

( )

(

)

]

(

k

)

k

( )

t k k t t k k k t _{k k k} k 1 1 1 2 1 3 2 2 2 1 1 2 3 2 − + = + ₊ + φ − φ + α − ₊ + ₋ φ φ (D4) PoniewaŜ

(

)

[

]

= −

(

+

) (

∏

[

+

)

−

]

∏ + − − − = − = 2 0 1 0 1 2 2 1 2 k i k i t i t t i α α α α i podobnie

(

)

[

]

₍

₎

∏

[

(

+

)

+

]

+ = ∏ + + + + = = 1 0 0 1 2 1 2 1 2 k i k i t i t t i α α α α

na podstawie wzoru (D3) dla k otrzymujemy

(

)

(

)

(

)

[

(

)

]

[

(

)

]

(

)

(

)

[

k t

]

[

(

k

)

t

]

( )

t t t i t i k t t k k i k i k φ α α α α α α α α φ + + − − + − = ∏ + + ∏ + − + + − = + + = − = 3 2 1 2 1 2 1 2 1 2 2 2 2 1 0 2 0 2 (D5)

Następnie, na podstawie wzorów (D4) i (D5), uwzględniając formułę (D3) dla

(

k−1

)

i k , po prostych przekształceniach algebraicznych, otrzymujemy wzór (D3) dla

(

k+1

)

.

Obecnie przejdziemy do wykazania słuszności wzoru (6); wykorzystamy wzór (D1). Na mocy wzoru (D3) dla dowolnego k ≥ 0

( )

(

)

[

(

)

]

(

)

[

]

(

)

(

)

(

2 1

)

2 1 2 1 2 2 1 2 1 2 1 2 2 0 0 1 0 + = + + = + + + =

∏

∏

= − = k k k i i k k i k i k α α α α α α φ (D6)

a stąd na podstawie wzorów (D3) oraz (D1) otrzymujemy dla k ≥0

( )

(

)

(

)

(

)

[

]

(

)

[

]

∏

∏

= − = + + − + + − + = _k i k i k t i t i k k t 0 1 0 1 2 1 2 1 2 2 1 2 2 α α α α ϕ

skąd po prostych przekształceniach natychmiast wynika wzór (6).

Wzór (8) wynika z (6) dla k =0. Aby udowodnić słuszność reguły rekurencyjnej (7) wystarczy zauwaŜyć, Ŝe na podstawie wzoru (D3) dla kaŜdego k ≥0 zachodzi zaleŜność rekurencyjna

( )

( )

_[

[

₍

(

)

₎

]

_]

_{, k ,,}K t k t k k k t t _k k 01 3 2 1 2 1 2 3 2 1 ₊+ −₊ = + + = + φ α_α φ

Stąd, wykorzystując wzór (D1), otrzymujemy równość

( )

( )

_[

[

₍

(

)

₎

_]

]

( )

_[

[

₍

(

)

₎

]

_]

t k t k k k t k k k k t _{k k} k ₊+ −₊ + + − + + + + = + φ φ α_α ϕ 3 2 1 2 1 2 3 2 3 2 1 2 1 2 3 2 0 1

czyli, wykorzystując ponownie wzór (D1), mamy

( )

( )

_[

[

₍

(

)

₎

]

_]

[

( )

( )

]

_[

[

₍

(

)

₎

]

_]

t k t k k k t k k k k t _{k k k} k ₊+ −₊ + + − − + + + + = + φ φ ϕ _αα ϕ 3 2 1 2 1 2 3 2 0 3 2 1 2 1 2 3 2 0 1

a stąd, poniewaŜ na podstawie wzoru (D3) φ_k

( )

0 = 2α 2k+1, otrzymujemy

( )

( )

_[

[

₍

(

)

₎

]

_]

(

₍

₎

)

_[

[

₍

(

)

₎

]

_]

_      + + − + + + − + + + + − + + + = + k_{k k}k t_t k t k t k k k t t _k k ϕ α_α α _αα ϕ 3 2 1 2 1 2 3 2 1 3 2 2 3 2 1 2 1 2 3 2 1

PIŚMIENNICTWO

1. Bubnicki Z.: Identification of Control Plants. PWN, Warszawa, Elsevier, Amsterdam, 1980. 2. Bzowska-Bakalarz M.: Właściwości mechaniczne korzeni buraków cukrowych. Rozprawy

Naukowe Akademii Rolniczej w Lublinie, 166, 1994.

3. Christensen R. M.: Theory of Viscoelasticity. An Introduction. Academic Press, New York, 1971. 4. Elster C., Honerkamp J., Weese J.: Using regularization methods for the determination of relaxation and retardation spectra of polymeric liquids. Rheological Acta, 30, 161-174, 1991. 5. Fujihara S., Yamamoto R., Masuda Y.: Maxwellian Spectra of Stress Relaxation in the Cell

Wall and Growth Regulation in Higher Plants. Proc. of Int. Workshop Stress Relaxation in Solids and Biological Origin, 47-51, Prague, 1995.

6. Golub G.H., Heath M.T., Wahba G.: Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21, 215-223, 1979.

7. Gu Ch., Heckman N., Wahba G.: A note on generalized cross-validation with replicates. Statistics & Probability Letters, 14, 283-287, 1992.

8. Hellebrand H. J.: Comparison of Models for Evaluation of Stress Relaxation. Proc. of Int. Workshop Stress Relaxation in Solids and Biological Origin, 3-10, Prague, 1995.

9. Kaczorek T.: Wektory i macierze w automatyce i elektrotechnice. WNT, Warszawa, 1998. 10. Lebiediew N.N.: Funkcje specjalne i ich zastosowania. PWN, Warszawa, 1957.

11. Lemaitre J. (Ed.): Handbook of Materials Behavior Models, Academic Press, New York, 2001. 12. Malkin A. Ya., Masalova I.: From dynamic modulus via different relaxation spectra to relaxation

and creep functions. Rheol. Acta, 40, 261-271, 2001.

13. Nguyen N., Milanfar P., Golub G.: A computationally efficient superresolution image reconstruction algorithm. IEEE Trans. Image Processing, 10(4), 573-583, 2001.

14. Owens R. G., Phillips T. N.: Computational Rheology. Imperial College Press London, 2002. 15. Paulson K. S., Jouravleva S., McLeod C. N.: Dielectric Relaxation Time Spectroscopy. IEEE

Trans. on Biomedical Engineering, 47, 1510-1517, 2000.

16. Ptaszek P., Grzesik M.: LepkospręŜyste zjawiska retardacyjne w modelowym układzie skrobia kukurydziana-guma guar. Acta Agrophysica, 3(3), 573-578, 2004.

17. Rao M.A.: Rheology of Fluid and Semisolid Foods. Principles and Applications. Aspen Publishers, Inc., Gaithersburg, Maryland, 1999.

18. Stankiewicz A.: Algorytm identyfikacji ciągłego spektrum czasów relaksacji biologicznych materiałów lepkospręŜystych. Acta Sci. Pol., Seria Technica Agraria, 2(2), 77-91, 2003. 19. Stankiewicz A.: Metoda optymalnej identyfikacji ciągłego spektrum relaksacji materiałów

lepkospręŜystych. InŜynieria Rolnicza, 12(54), 373-384, 2003.

20. Stankiewicz A., Gołacki K.: Metoda wyznaczania ciągłego spektrum relaksacji lepkospręŜystych materiałów roślinnych. Acta Agrophysica, 2(3), 627-638, 2003.

21. Ter Haar D.: A Phenomenological Theory of Visco-Elastic Behaviour. I. Physica, XVI, 719-737, 1950.

22. Tikhonov A. N., Arsenin V. Y.: Solutions of Ill-posed Problems. Wiley, New York, 1977. 23. Wahba G., Wold S.: A completely automatic French curve: fitting splines by cross validation.

APPROXIMATION OF THE RETARDATION SPECTRUM

OF VISCOELASTIC PLANT MATERIALS USING LAGENDRE FUNCTIONS

Anna Stankiewicz

Department of Technical Science, University of Agriculture ul. Doświadczalna 50A, 20-280 Lublin

e-mail: anna.stankiewicz@ar.lublin.pl

A b s t r a c t . The paper deals with the problem of direct recovery of continuous retardation spectrum of linear viscoelastic materials from discrete-time noise measurements of creep compliance obtained in a standard creep experiment. An optimal orthogonal scheme of the spectrum approximation by the finite series of Legendre functions is presented. Since the problem of retardation spectrum identification is practically an ill-posed problem of reconstructing solution of Fredholm integral equation of the first kind from measured data, Tikhonov regularization with generalized cross validation (GCV) is used to guarantee the stability of the scheme. It is proved that the accuracy of the spectrum approximation depends both on measurement noises and regularization parameter as well as on the proper selection of the basic orthogonal functions. Numerical calculations on model data are enclosed.

K e y w o r d s : viscoelasticity, retardation spectrum, identification algorithm, regularization, Legendre functions