Tomasz HANISZEWSKI[1]

METODYKA IDENTYFIKACJI PARAMETRÓW ZMODYFIKOWANEGO MODELU BOUCA-WENA

Streszczenie. Model Bouca-Wena jest sformułowaniem teoretycznym, umożliwiającym odzwierciedlenie rzeczywistej pętli histerezy modelowanego obiektu, jakim może być lina nośna występująca w mechanizmie podnoszenia suwnicy, gdzie dla przyjętej, zmodyfikowanej wersji modelu występuje dziewięć parametrów. Wyznaczenie takiej liczby parametrów jest zagadnieniem złożonym i problematycznym. W niniejszym artykule przedstawiono zatem metodykę identyfikacji oraz przykładowe wyniki symulacji numerycznych. Otrzymane wyniki porównano z danymi uzyskanymi na podstawie badań laboratoryjnych lin [3] i na ich podstawie stwierdzono zgodność wyników i możliwość zastosowania w układach dynamicznych zawierających w swej strukturze cięgna [4].

Słowa kluczowe: Model, Bouc-Wen, lina, modelowanie, identyfikacja

methodology for the estimation of parameters, of the modified Bouc-Wen model

Summary. Bouc-Wen model is theoretical formulation that allows to reflect real hysteresis loop of modeled object. Such object is for example a wire rope, which is present on equipment of crane lifting mechanism. Where adopted modified version of the model has nine parameters. Determination of such a number of parameters is complex and problematic issue. In this article are shown the methodology to identify and sample results of numerical simulations. The results were compared with data obtained on the basis of laboratory tests of ropes [3] and on their basis it was found that there is compliance between results and there is possibility to apply in dynamic systems containing in their structures wire ropes [4].

Keywords: Model, Bouc-Wen, wire rope, modeling, identification

1. WPROWADZENIE

Mianem histerezy w naukach przyrodniczych określane jest zjawisko zależności aktualnego stanu układu od stanów w poprzedzających chwilach [2]. Jego graficznym przedstawieniem jest pętla, a samo zjawisko jest nieliniowe i występuje w wielu systemach fizycznych. Mogą to być systemy z tarciem, jak na przykład liny będące jednym
z podstawowych elementów maszyn transportowych, takich jak np. dźwignice.

Histereza występuje w układach z tarciem, a układ taki przedstawia między innymi charakterystyka siłowo-przemieszczeniowa, uzyskana na podstawie badań lin, co opisano
szczegółowo w pozycjach [3, 4, 10], gdzie z wykresów rozciągania wynika, że zależność pomiędzy obciążeniem a wydłużeniem jest krzywoliniowa, a przebieg krzywej podczas odciążania próbki jest inny niż podczas obciążania.

W celu modelowania zjawiska histerezy występującego w cięgnie, należy zastosować konkretny model matematyczny, który określa właściwości mechaniczne liny w postaci odwzorowania pętli histerezy. Wiąże się z tym określenie występującego w linie nieliniowego zjawiska dyssypacji energii oraz nieliniowej sztywności, co stwarza pewne trudności. Prezentowany w artykule model fenomenologiczny pozwala na wprowadzenie właściwości mechanicznych liny na potrzeby symulacji numerycznej, np. procesu podnoszenia ładunku
w suwnicach, i jest alternatywą metod modelowania opartych na metodzie regresji
i równań wielomianowych co najmniej trzeciego rzędu, lub do dotychczas stosowanych modeli, np. dwójnika sprężysto-tłumiącego Kelvina-Voigta. Zatem głównym celem wprowadzenia modelu Bouca-Wena jest zwiększenia dokładność obliczeń, zwłaszcza podczas wyznaczenia obciążeń dynamicznych w mechanizmach podnoszenia suwnic, poprzez uwzględnienie w modelu numerycznym liny zjawiska histerezy.

2. MODEL FENOMENOLOGICZNY BOUCA-WENA

W literaturze [10, 11, 12, 13, 14, 15, 18] opisano liczną grupę modeli matematycznych, umożliwiających przeprowadzenie numerycznej symulacji histerezy. Modele te różnią się głównie stopniem skomplikowania, a dokładniej liczbą estymowanych parametrów oraz złożonością matematyczną. Wśród nich można wymienić między innymi takie modele, jak: modele Duhema, Preisacha, Bouca-Wena czy Coulomba. W celu możliwie dokładnego odwzorowania zjawiska histerezy, zachodzącego w linie nośnej badanego obiektu wybrano model zmodyfikowany Bouca-Wena, gdyż jest on bardzo szeroko stosowany do opisu zjawiska histerezy i w swej strukturze zawiera nieliniową sztywność i odpowiednią funkcję, odpowiedzialną za możliwość kształtowania pętli. Proponowany model cechuje się dużą swobodą w opisie układów z histerezą, a zarazem od strony matematycznej jest on relatywnie prosty. Za jego wadę można uznać liczbę parametrów potrzebnych do opisu modelu. Umożliwia on zapis w postaci analitycznej i przedstawia zachowanie szerokiej klasy układów, gdzie występuje zjawisko histerezy. Jednym z głównych problemów przy wykorzystaniu omawianego modelu, jest (jak wspomniano) problematyka identyfikacji parametrów. Model ten przedstawia się za pomocą następującej zależności [6, 7, 8, 13]:

(1)

W artykule zdecydowano się na użycie modelu zmodyfikowanego, tj. rozbudowanego
o połączenie tzw. siły zwrotnej opartej na modelu podstawowym Bouca-Wena z układem nieliniowej sprężystości i funkcji modulującej kształt pętli. Zmodyfikowany model
Bouca-Wena przedstawiono poniżej [13]:

(2)

Jak opisuje układ równań, całkowita siła generowana przez model wynika z połączenia równania, opisującego nieliniową sztywność (funkcja F₁(t)) oraz funkcji wykładniczej F₂(t), będącej formą funkcji modyfikującej kształt pętli. Wszystkie 9 parametrów modelu wraz
z jednostkami przedstawiono w tablicy 1. Parametry modelu zestawiono w postaci wektora (3), który w dalszej części posłuży do symulacji i optymalizacji.

(3)

Tablica 1 Parametry przyjętego zmodyfikowanego modelu Bouca-Wena
Lp.	Parametr	Wymiar	Lp.	Parametr	Wymiar
1	α	[Nm^-1]	6	c	[m^-1]
2	β	[N^(1-n)m^-1]	7	k1	[Nm^-1]
3	γ	[N^(1-n)m^-1]	8	k2	[Nm^-2]
4	n	[–]	9	k3	[Nm^-3]
5	b	[–]

Zmienne k₁, k₂, k₃ określają nieliniową funkcję sztywności, a ponieważ funkcja F₂(t) jest bezwymiarową funkcją modulującą kształt, więc wymiary współczynników b oraz c będą bezpośrednio związane z równaniem określającym siłę F₂(t). Dokładny opis wpływu parametrów na zachowanie modelu oraz na genezę matematyczną można znaleźć
w literaturze [6, 9, 12, 13, 17, 18]. Układ prezentowanych równań (2) zamodelowano
w postaci schematów blokowych Matlab-Simulink, które posłużą jako modele do przeprowadzenia symulacji numerycznej.

3. METODYKA EKSPERYMENTU NUMERYCZNEGO

Równania fenomenologicznego modelu Bouca-Wena zawierają dziewięć parametrów determinujących kształt pożądanej charakterystyki. Problem identyfikacji można sprowadzić do problemu optymalizacji i za pomocą eksperymentu numerycznego znaleźć poszukiwane wartości. Gdy mamy wyniki badań doświadczalnych liny oraz możliwość przeprowadzenia symulacji numerycznej rozpatrywanego obiektu, odwzorowując jego zachowanie, można stwierdzić, że na podstawie porównania krzywych eksperymentalnych oraz numerycznych można wyznaczyć wartości parametrów modelu, kontrolując różnicę pomiędzy tymi krzywymi w zależności od parametrów modelu, gdzie jako funkcję celu przyjęto sumę błędów w określonych punktach czasowych badanego zakresu. Rozwiązując zatem zadanie minimalizacji różnic, można wyznaczyć wartości parametrów modelu. Problem optymalizacji w przypadku badanego systemu zmodyfikowany model Bouca-Wena można zapisać
w postaci problemu minimalizacji zadanej funkcji celu:

, (4)

gdzie:

Y – wektor wartości ze zbioru wartości dopuszczalnych,

– zbiór wartości dopuszczalnych,

– wartość siły w danym punkcie uzyskana za pomocą symulacji,

– wartość siły w danym punkcie uzyskana za pomocą badań doświadczalnych,

n – liczba iteracji.

Celem minimalizacji jest znalezienie wektora parametrów Y_min, minimalizującego funkcje celu. Rozpatrywany problem optymalizacji można także zapisać w postaci równania (5), biorąc pod uwagę układ równań (2). Zależność (5) powstała z przekształcenia układu równań (2). Otrzymana zależność w postaci równania (5) będzie zawsze spełniona, tj. będzie równa zeru, zatem, wystarczy zamienić siłę generowaną przez model F(t) siłą otrzymaną na drodze badań doświadczalnych, tj. F_bad (6). W momencie gdy parametry modelu będą zbliżać się do optymalnych wartość, funkcja określona równaniem (6) będzie zmierzała do zera. Dlatego też, optymalizując parametry modelu dąży się do minimalizacji wartości funkcji (6).

(5)

(6)

W rozważanym zagadnieniu różnice pomiędzy krzywymi eksperymentalną a numeryczną wyrażono, porównując wartości sił zarejestrowanych doświadczalnie i sił obliczonych za pomocą modelu Bouca-Wena, w punktach czasowych odpowiadających jednemu okresowi wymuszenia. W rozważanym zagadnieniu szukany jest zatem wektor parametrów Y minimalizujący funkcję (4) lub (6).

Do rozwiązania problemu wybrano algorytm Levenberga-Marquardta (LMA), który jest jednym z częściej wykorzystywanych algorytmów optymalizacji nieliniowej, a głównie
w nieliniowym zadaniu najmniejszych kwadratów. Jako algorytm iteracyjny łączy cechy metody największego spadku i Gaussa-Newtona. Dodatkowo, wykorzystano metodę Pattern Search (PS), czyli wyszukiwanie wzorca, należącą do grupy metod numerycznych optymalizacji. Podejście to wynikało z obserwacji, ponieważ podczas symulacji wstępnych zauważono, że początkowe wykorzystanie metody (PS), czyli wstępne dopasowanie funkcji, powoduje, że algorytm (LMA) na podstawie metody nieliniowej najmniejszych kwadratów uzyskuje zbieżność przy mniejszej liczbie iteracji. Dlatego też zdecydowano się na estymację dwuetapową, najpierw z użyciem metody (PS), a następnie na podstawie uzyskanych wartości parametrów ich optymalizację za pomocą algorytmu (LMA). W przypadku (PS) oraz (LMA) wskaźnikiem błędu była wartość tzw. funkcji błędu, co w oprogramowaniu Maltab-Simulink nazywane jest funkcją kosztu SSE [1, 16] lub celu, będącą niczym innym, jak sumą kwadratów odchyłek, co opisano we wzorze funkcji (4).

Do przeprowadzenia eksperymentu numerycznej estymacji parametrów badanego modelu wykorzystano narzędzie (SDO) „Simulink Design Optimization”, będące składnikiem pakietu Matlab-Simulink. „Simulink Design Optimization” umożliwia zwiększenie dokładności odwzorowania modelu dzięki porównywaniu z danymi testowymi.

W celu przeprowadzenia eksperymentu numerycznego, zbudowano model zawierający system wejścia i wyjścia dla narzędzia (SDO). Schemat blokowy przedstawiono na rys. 1. Głównym problemem przy wykorzystaniu (SDO) było odczytywanie danych z badań, ponieważ były one zapisywane jedynie w postaci wydruku. Problem ten jest bezpośrednio związany z brakiem przystosowania stanowiska badawczego do współczesnych systemów obliczeniowych, gdzie niezbędna jest możliwości eksportu danych w postaci ciągu wartości liczbowych. W związku z tym, została napisana aplikacja w postaci skryptu pakietu Matlab, umożliwiająca odczyt punktów pomiarowych na bazie analizy numerycznej obrazu,
a problem ten opisano dokładnie w pracy [5]. Zastosowanie programu do konwersji obrazów pozwoliło na uzyskanie chmury punktów pomiarowych do dalszej obróbki danych. Wykorzystanie opisanej metody analizy obrazu pozwoliło skrócić czas potrzebny do przygotowania danych do edycji, przy zachowaniu dostatecznej dokładności.

Wykorzystując przetworzone dane pomiarowe, przeprowadzono dwuetapową operację estymacji i optymalizacji parametrów. Proces ten polegał na imporcie danych wejściowych
w postaci wymuszenia użytego podczas badań laboratoryjnych lin w dziedzinie czasu oraz imporcie danych wyjściowych uzyskanych podczas badań lin, w tym przypadku siły
w dziedzinie czasu. W modelu przyjęto, że wartości początkowe dla każdego parametru wynoszą 0. Tak przygotowany system umożliwia estymację, a kolejnym etapem jest ustawienie kroku całkowania oraz rodzaju solvera. W rozpatrywanym przypadku wybrano solver stałokrokowy ode5 (Dormand-Prince).

Rys. 1. Schemat blokowy zmodyfikowanego modelu Bouca-Wena

Fig. 1. Block diagram of the modified Bouc-Wen model

Warto podkreślić, że określenie najbardziej optymalnego solvera dla konkretnego modelu na ogół wymaga eksperymentowania. Przy wyborze solvera stałokrokowego w dokumentacji podaje się pewną metodykę eksperymentowania, proponując wybór kolejno solvera ode1, ode3, ode2, ode5, ode4 [1, 16]. Po dokonaniu wyboru solvera wstępnie wyselekcjonowano metodę (PS) i przeprowadzono estymację, po czym użyto wyznaczonych parametrów
tzw. inicjujących w celu ich optymalizacji metodą nieliniową najmniejszych kwadratów,
z wykorzystaniem algorytmu (LMA). Przy wyznaczaniu serii parametrów dotyczących różnych danych związanych z cyklem obciążenia i wartością obciążeń, metodę (PS) należy zastosować tylko podczas wyznaczania pierwszego zestawu parametrów, tzw. inicjujących. Następnie zestawy estymowanych danych można wykorzystać do optymalizacji tylko przy użyciu algorytmu (LMA).

4. WYNIKI EKSPERYMENTU NUMERYCZNEGO

Na podstawie przeprowadzonych eksperymentów numerycznej estymacji i optymaliacji parametrów określono wartości niezbędnych parametrów modelu. Czas potrzebny na estymację oraz uzyskanie pozostałych informacji dla kolejnych partii danych przedstawiono w tablicy 2, gdzie zestawiono przykładowe wyniki.

Tablica 2 Porównanie przebiegów estymacji dla przyjętego zmodyfikowanego modelu Bouca-Wena
A	B	C	D	E	F	G	H
m=6.3	1	82	36	232	1500	0,05	PS
m=6.3	1	662	226	227	300	0,05	LMA
m=6.3	1	82	36	232	1500	0,05	PS
m=6.3	10	1856	399	400	550	0,01	LMA
m=2	10	1365	156	158	450	0,01	LMA
m=1.5	10	3360	339	340	920	0,01	LMA

gdzie: A – współczynnik bezpieczeństwa, B – cykl, C – czas [s], D – liczba iteracji, E – liczba wywołań funkcji, F – maksymalna różnica względem sygnału referencyjnego [N], G – krok całkowania, H – algorytm, metoda (Levenberga-Marquardta (LMA), Pattern search” (PS)).

W tablicy 3 przedstawiono przykładowe, wyznaczone wartości parametrów dla wybranych cykli w zależności od współczynnika bezpieczeństwa i cyklu obciążenia, gdzie współczynnik bezpieczeństwa określony jest jako stosunek rzeczywistej siły zrywającej linę do siły odpowiadającej danemu obciążeniu.

Tablica 3 Wartości estymowanych parametrów przyjętego zmodyfikowanego modelu Bouca-Wena
Lp.	Parametr	Cykl 01 dla m=6.3	Cykl 10 dla m=6.3	Cykl 01 dla m=2	Cykl 10 dla m=2
1	k1	-3,77E+03	-1,22E+05	-9,26E+04	-1,04E+05
2	k2	4,87E+04	4,84E+04	4,71E+04	4,72E+04
3	k3	2,15E+09	2,15E+09	2,15E+09	2,15E+09
4	b	1,93E+10	1,93E+10	1,93E+10	1,93E+10
5	c	-7,35E+00	-7,53E+00	-3,27E+00	-2,35E+00
6	α	2,88E+05	1,92E+05	3,58E+05	2,13E+05
7	β	3,46E+02	2,50E+02	3,11E+02	6,21E+02
8	γ	-2,89E+03	-3,22E+03	-4,55E+03	-1,23E+04
9	n	8,20E-01	8,13E-01	7,18E-01	6,05E-01

Dla otrzymanych wartości parametrów zasymulowano układ z rys. 1 i przedstawiono porównanie danych otrzymanych na postawie symulacji oraz analizy obrazu i odczytu na maszynie wytrzymałościowej. Przykładowo, na rys. 2, 3 przedstawiono wykresy ukazujące porównanie danych dla pierwszego oraz dziesiątego cyklów obciążenia, przy współczynniku bezpieczeństwa 2 oraz 6.3, co odpowiada sile: 45 [kN] i 15 [kN] maksymalnego obciążenia
w stosunku do wartości siły powodującej jej zerwanie.

a) b)

∙10^-3

∙10^-3

Rys. 2. Wykres ujmujący zależność przemieszczenia od siły dla danych, uzyskanych na podstawie analizy obrazu, odczytu wskazań maszyny oraz symulacji dla m=2 (a – cykl 1, b – cykl 10)

Fig. 2. Graph showing the relationship between the force and displacement for data obtained by image analysis, the reading and the simulation for m=2 (a – cykl 1, b – cykl 10)

a) b)

∙10^-3

Rys. 3. Wykres ujmujący zależność przemieszczenia od siły dla danych, uzyskanych na podstawie analizy obrazu, odczytu wskazań maszyny oraz symulacji dla m=6.3 (a – cykl 1, b – cykl 10)

Fig. 3. Graph showing the relationship between the force and displacement for data obtained by image analysis, the reading and the simulation for m=6.3 (a – cykl 1, b – cykl 10)

5. WNIOSKI

Jak można zaobserwować na rys. 2, 3, przebiegi reprezentowane przez symulacje za pomocą układu równań Bouca-Wena z użyciem wyznaczonych wartości parametrów są bardzo zbliżone do przebiegów wyznaczonych podczas badań, co świadczy o dobrym dopasowaniu danych do modelu.

Model zmodyfikowany Bouca-Wena, zastosowany do opisu cięgna umożliwia zatem wprowadzenie właściwości liny do modelu dynamicznego, usprawniając badania modelowe, a jednocześnie ułatwiając i upraszczając proces analizy oraz wyznaczania jej parametrów fizycznych. Zatem, proponowany model można zastosować w układach służących symulacji procesu podnoszenia ładunku, np. w suwnicach, a dokładnie odwzorowana charakterystyka obciążeniowo-odciążeniowa, reprezentowana przez model Bouca-Wena, umożliwi badanie wpływu zastosowanego modelu cięgna na wartości współczynnika nadwyżki dynamicznej,
w odniesieniu np. do wartości normatywnych lub otrzymanych na podstawie zastosowania klasycznego dwójnika sprężysto-tłumiącego.

Badania były współfinansowane ze środków Unii Europejskiej w ramach Europejskiego Funduszu Społecznego, w projekcie „Aktywizacja społeczności akademickiej jako element realizacji Regionalnej Strategii Innowacji” POKL.08.02.01-24-019/08”. Praca była wykonana jako cześć BKM-297/RT3/2012 i BKM-511/RT3/2013.

Bibliografia

1. Brzózka J., L. Dorobczyński. 2008. Matlab – środowisko obliczeń numerycznych. [In Polish: Matlab - Environment numerical calculations]. Warszawa: Mikom - PWN.

2. Encyklopedia popularna PWN. [In Polish: Encyclopedia PWN]. 2008. Wydanie trzydzieste czwarte, zmienione i uzupełnione. Warszawa.

3. Haniszewski T., D. Gąska. 2011. „Line 6x19 Seale +fc sZ hysteresis determination”. Zeszyty Naukowe Politechniki Śląskiej, seria Transport 73: 21-30. Gliwice.

4. Haniszewski T. 2013. „Modelowanie dynamiki lin stalowych w konstrukcjach maszyn transportowych”. Rozprawa doktorska. [In Polish: “Modeling dynamic ropes in the construction machinery transport”. PhD thesis]. Katowice: Politechnika Śląska.

5. Haniszewski T. 2011. „Conversion of experimental charts to digital data as a tool for further processing of results”. Dziennik Wschodnioukraińskiego Narodowego Uniwersytetu 12(166): 250-257.

6. Ikhouane F., V. Mañosa, J. Rodellar. 2007. “Dynamic properties of the hysteretic Bouc-Wen model”. Systems and Control Letters 56(3): 197-205.

7. Ikhouane F., J. Rodellar. 2007. Systems with Hysteresis - Analysis, Identification and Control using the Bouc–Wen Model. UK, Chichester.

8. Ismail M., F. Ikhouane, J. Rodellar. 2009. “The hysteresis bouc-wen model, a survey”. Archives of Computational Methods in Engineering 16(2): 161-188.

9. Krasnosel'skii M.A., A.V. Pokrovskii. 1989. Systems with Hysteresis. New York: Springer-Verlag.

10. Margielewicz J., T. Haniszewski, D. Gąska, C. Pypno. 2013. Badania modelowe mechanizmów podnoszenia suwnic. Wydawnictwo J&L Leszek Żochowski. Katowice.

11. Mayergoyz I.D. 1991. Mathematical Models of Hysteresis. New York: Springer-Verlag.

12. Ni Y.Q., J.M. Ko, C.W. Wong, S. Zfaan. 1999. “Modelling and identification of a wire-cable vibration isolator via a cyclic loading test part 1: experiments and model development”. Proceedings of the Institution of Mechanical Engineers, part I: “Journal of Systems and Control Engineering” 213(3): 163-171.

13. Nijmeijer H., R.H.B. Fey, J.E. van Aanhold, M.A. Hulsen, H. Nijmeijer. 2004. Modelling and identification of dynamic behaviour of a wire rope spring. Technische Universiteit Eindhoven, The Netherlands, Department of Mechanical Engineering, Section Dynamics and control technology, Eindhoven, 1st March 2004. Report No. DCT-2004/28.

14. Oh J., A.K. Padthe, D.S. Bernstein, D.D. Rizos, S.D. Fassois. 2005. “Duhem models for hysteresis in sliding and presliding friction. In Proceedings of the 44th IEEE Conference on Decision and Control, and the European Control Conference, CDC-ECC '05 2005, No. 1583478. P. 8132-8137.

15. Preisach F. 1935. „Über die magnetischenachwirkung“. [In German: „Use the magnetic aftereffect“]. Zeitschriftfür Physik 94(5-6): 277-302.

16. Skowronek M. 2008. Modelowanie cyfrowe. [In Polish: Digital modeling]. Gliwice: Wydawnictwo Politechniki Śląskiej.

17. Thyagarajan R.S. 1989. Modeling and analysis of hysteretic structural behavior. Earthquake engineering research laboratory, CaltechEERL. Report No. EERL-89-03.

18. Visintin A. 1994. Differential Models of Hysteresis. Berlin: Springer.

---

[1] Silesian University of Technology, Faculty of Transport, Katowice, Poland, e-mail: tomasz.haniszewski@polsl.pl