Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 UDK - UDC 517.518.8:519.243 Strokovni članek - Speciality paper (1.04) Uporaba metode premičnih najmanjših kvadratov za gladko aproksimacijo vzorčenih podatkov The Use of Moving Least Squares for a Smooth Approximation of Sampled Data Igor Grešovnik (Center za računalništvo v mehaniki kontinuuma, Ljubljana) Predstavljena je uporaba metode premičnih najmanjših kvadratov za glajenje podatkov in aproksimacijo odzivnih funkcij s šumom. Obravnavane so lastnosti aproksimacije glede na raven šuma, gostoto vzorčenja in dejanski doseg vpliva vzorcev. Uporabnost metode je prikazana pri glajenju merskih podatkov in reševanju optimizacijskega problema s šumom z metodo zaporednih aproksimacij. Metoda premičnih najmanjših kvadratov se izkaže za uporabno pri reševanju različnih problemov zaradi gladkosti, možnosti natančne aproksimacije na poljubno velikem območju, možnosti aproksimacije na podlagi neurejene množice vzorčnih točk in prilagodljivosti različnim lastnostim aproksimirane funkcije v različnih območjih. © 2007 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: vzorčni podatki, aproksimacije funkcij, metode premičnih najmanjših kvadratov, glajenje, optimiranje) The use of the moving least-squares approximation for the smoothing of data and the approximation of noisy response functions is presented. The approximation properties are treated with respect to the level of noise, the sampling density and the effective range of the sample. The applicability of the method is demonstrated by smoothing experimental measurement data and solving an optimization problem with a noisy response by the successive response approximation technique. The results indicate that the moving least-squares approximation is applicable to a wide variety of problems due to its smoothness, its accurate approximation over arbitrarily large domains using low-order basis functions, its ability to deal with an irregular arrangement of sampling points and to adapt to different modes of the approximation function in different regions. © 2007 Journal of Mechanical Engineering. All rights reserved. (Keywords: sample data, function approximation, moving least squares methods, smoothing, optimization) 0 UVOD 0 INTRODUCTION Pogosto želimo množico točk, ki predstavljajo neko funkcijsko zvezo, zamenjati s približno funkcijo z določenimi zveznostnimi lastnostmi. Ena od možnosti je aproksimacija s polinomi določenega reda ali s trigonometrijsko vrsto. To je učinkovito, kadar potrebujemo približek na omejenem območju. V nasprotnem primeru potrebujemo veliko število členov aproksimacije, posebej, ko imamo več neodvisnih spremenljivk. Aproksimacija s polinomi je nestabilna, ko je število točk veliko in se število členov približa številu točk. Ta problem lahko rešujemo z odsekovno polinomsko aproksimacijo [4]. Pri tem postopku se odpovemo zveznosti poljubne stopnje. V primeru več neodvisnih spremenljivk navadno It is often useful to replace a set of points that represent a functional relation with an approximate relation that has certain continuity properties. One of the possibilities is an approximation with a polynomial of a given order or a series of trigonometric functions. This approach is efficient when an approximation is needed over a limited domain. In other cases a large number of terms is necessary, especially in the multivariate case. The polynomial approximation also becomes unstable when the number of points is large and close to the number of terms. However, this problem can be tackled by a piecewise polynomial approximation [4]. With this approach, one gives up the continuity of an arbitrary order and in the 582 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 potrebujemo strukturirano razdelitev območja aproksimacije. Kot druga možnost je v tem prispevku predstavljena metoda premičnih najmanjših kvadratov. Metoda omogoča gladko in stabilno aproksimacijo na poljubno velikem območju. Ni potrebna kakršnakoli delitev območja aproksimacije in za gladke funkcije lahko dosežemo poljubno natančnost aproksimacije z uporabo omejenega števila osnovnih funkcij, če lahko ustrezno povečujemo gostoto vzorčenja. Oris metode je podan v poglavju 1. Lastnosti aproksimacije so v poglavju 2. prikazane na primeru z analitično funkcijo ene spremenljivke. Dva primera praktične uporabe sta prikazana v poglavju 3. V prvem primeru je metoda uporabljena za glajenje časovno odvisnega signala zaznavala pri laboratorijski meritvi, kar omogoči učinkovito uporabo minimizacijskega postopka pri določitvi modelskih parametrov z reševanjem obrnjenega problema. V drugem primeru je metoda uporabljena pri zaporednih prilagodljivih aproksimacijah namenske in omejitvenih funkcij. Aproksimacije uporabimo v iterativnem optimizacijskem postopku, primernem za reševanje problemov s šumom. multivariate case a structured division of the domain of approximation is usually needed. As an alternative, the moving least-squares approximation method is presented. The method enables a smooth and stable approximation over arbitrarily large domains. No particular partition of the domain of approximation is necessary and arbitrary accuracy can be achieved for smooth functions with a small number of basis functions, provided that the sampling density can be increased correspondingly. An outline of the method is given in Section 1. The properties of the approximation are studied in Section 2 on a synthetic example with a function of a single variable. In Section 3, two practical applications of the method are demonstrated. In the first example the method is used for the smoothing of a time-dependent sensor output in a laboratory experiment, which makes it possible to perform an efficient minimization procedure for an inverse estimation of the model parameters. In the second example the method is used for successive adaptive approximations of the noisy objective and constraint functions, which are used in an iterative optimization procedure suitable for optimization with a noisy response. 1 OPIS METODE 1 METHOD DESCRIPTION 1.1 Linearna aproksimacija po metodi najmanjših kvadratov z utežmi 1.1 Linear Weighted Least-Squares Approximation Pri linearni metodi najmanjših kvadratov aproksimiramo neznano funkcijo f{x) z linearno kombinacijo n osnovnih funkcij ^(x), ..., fjx) na podlagi znanih (vzorčenih) vrednosti funkcije v določenem številu točk: Using the least-squares method we approximate an unknown function f(x) by a linear combination of n basis functions f1(x), ..., fn(x) on the basis of known (sampled) values of the function at a number of points: (1). yk=f(xk) + rk, k = \,...,m Člen r. pomeni naključno napako (šum) pri merjenju ali računanju vrednosti funkcije. Aproksimacija: The term ri accounts for a random error (noise) that is eventually accomplished when the function is measured or evaluated. We want the approximation: j=1 (2) se mora čim bolje ujemati z vzorčenimi vrednostmi, torej: to agree as much as possible with the sampled values, i.e.: y{\)~yk Vk=u, m (3). V enačbi (2) so a. neznani koeficienti aproksimacije, ki jih je treba določiti. To storimo In Equation (2), aj are the coefficients of the approximation that must be determined. This is n Uporaba metode premičnih najmanjših kvadratov - The Use of Moving Least Squares 583 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 tako, da poiščemo najboljše ujemanje v skladu najmanjših kvadratov, torej z minimizacijo naslednje funkcije koeficientov: done by looking for the best agreement in the least-squares sense, i.e., by a minimization of the following function of coefficients: ±<{y{xk)-ykY = ± yk (4). V zgornji enačbi smo koeficiente a. uredili v vektor a. Nenegativne uteži wk merijo sorazmerno pomembnost vzorčnih točk. Večje so vrednosti uteži, bolj je aproksimacija prilagojena pripadajočim vzorčenim vrednostim na račun slabšega ujemanja z vrednostmi z manjšimi utežmi. Metoda najmanjših kvadratov z utežmi ima statističen pomen [1]. Če so merske napake rt porazdeljene normalno z znanimi standardnimi odstopanji a. in je (2) pravilen model za/(x), dobimo z minimizacijo 0(a) vrednosti koeficientov a z največjo verjetnostjo, da so izračunani koeficienti pravilni, če postavimo uteži na: Največja verjetnost se tu nanaša na vrh verjetnostne gostote za porazdelitev koeficientov a, ki jih dobimo z minimizacijo 0(a) pri različnih udejanjenjih izmerkov glede na njihovo verjetnostno porazdelitev. Čeprav porazdelitve merskih napak pogosto niso normalne in uporabljeni modeli niso natančni, je uporaba najmanjših kvadratov pri prilagajanju podatkov povsod navzoča in se izkaže za primerno v veliko primerih, ko nimamo na voljo fizikalno utemeljenih modelov. Minimizacijo 0(a) lahko izvedemo z iskanjem ustaljene točke, zato zahtevamo: Hajfj ( x In the above equation the coefficients a. were arranged in a vector a, and non-negative weighting coefficients wt measure the relative significance of the samples. The higher these coefficients are, the more the approximation will attempt to accommodate the corresponding sampled values on account of the poorer agreement with the samples with smaller weights. The weighted least-squares formulation has a statistical meaning [1]. When measurement errors r. have a normal distribution with known standard deviations G and (2) represents a correct model for /(x), a minimization of 0(a) yields the values of the coefficients a with the highest probability that they are correct, provided that the weights are set to: 1 °k (5). d(j)(a) da 21 The highest probability refers to the maximum of the probability density function for the distribution of coefficients a obtained for different realizations of measurements according to their probability distribution. Although the distributions of measurement errors are often not normal and in particular the model used is not correct, the least-squares approach is ubiquitously used for fitting data and proves suitable in many situations when we do not have physically founded models for the measured data at our disposal. The minimization of 0(a) is performed by finding the stationary point, i.e., by setting 7=1 )-yk /,( x*) 0 Vz = 1, (6). Dobimo sistem enačb za koeficiente a: This yields the system of equations for coefficients a: kjer so: in Ca = d where: k=1 and k=1 (7), (8) (9). 2 n k=1 w k 584 Grešovnik I. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 1.2 Premični najmanjši kvadrati 1.2 Moving Least Squares (MLS) Izbira osnovnih funkcij je odločilna za to, kako natančno lahko aproksimacijo prilagodimo podatkom. Če nimamo primernega fizikalnega modela, pogosto vzamemo množico enočlenikov do določene stopnje, kar je utemeljeno s Taylorjevim izrekom [2]. Za dobro aproksimacijo funkcije na večjem območju je treba ustrezno povečati število členov. Posebej pri večjem številu neodvisnih spremenljivk je lahko težko določiti primerno število osnovnih funkcij. Nekateri členi aproksimacije so lahko odveč kljub velikemu presežku podanih vrednosti glede na število koeficientov, zato je lahko sistem (7) slabo pogojen. Opisane težave lahko olajšamo s prostorsko omejitvijo aproksimacije, tako da je odvisna le od vzorčnih točk v omejeni okolici točke izračuna aproksimacije. Da ohranimo zveznost in zagotovimo zmožnost prilagoditve na celotnem območju, ki nas zanima, obdržimo obliko (2), vendar s spremenljivimi koeficienti, ki so zvezno odvisni od neodvisnih spremenljivk. Aproksimacijo lahko torej zapišemo kot: The choice of the basis functions has a crucial impact on the degree to which the approximation can accommodate the data. With the lack of an appropriate physical model, a set of monomials up to a given degree is often taken, with the justification relying on Taylor’s theorem [2]. When the function is to be approximated over a larger range of independent variables, an increased number of terms must be taken for a good approximation. In particular when the number of independent variables is large, it may be difficult to estimate a suitable number of basis functions. Some terms of the approximation may be redundant, and in spite of the large excess of data with respect to the number of coefficients, the system of Equations (7) can be badly conditioned. The above-mentioned difficulties can be alleviated by localizing the approximation in the sense that it depends only on samples in a certain neighborhood of the evaluation point. In order to preserve the continuity of the approximation and ensure its accommodation ability over the whole range of interest, we keep the form of the approximation similar to (2) but let coefficients continuously depend on independent variables. Therefore, we can write the approximation as: y(x) Koeficiente a(x) moramo posebej izračunati v vsaki točki, kjer izračunamo vrednost aproksimacije. V metodi premičnih najmanjših kvadratov [3] to dosežemo z uvedbo utežnih funkcij, ki se zmanjšujejo z razdaljo od pripadajočih vzorčnih točk. Koeficiente a(x) izračunamo posebej v vsaki točki računanja aproksimacije x z reševanjem sistema enačb, ki ustreza sistemu (7) do (9), pri čemer so uteži wt(x) odvisne od točke izračuna: Ya.(x)f.(x) (10)' 7=1 The coefficients a(x) are thus calculated at each point where the approximation is evaluated. In the moving least-squares method [3] this is done by introducing weight functions that fall with increasing distance from the corresponding samples. The system of equations equivalent to that defined by Equations (7) to (9) is then solved to obtain the coefficients a(x) in each evaluation point x, with weights wt(x) dependent on the point of evaluation: C(x)a(x) = d(x) C,(x) = I>(x)^(x,)/,(xJ (11). ^x=I^x!i(x) yk Utežne funkcije wk(x) se morajo v vseh smereh enakomerno zmanjševati z razdaljo od pripadajočih vzorčnih točk xt in morajo imeti določeno stopnjo zveznosti, s čimer zagotovimo zveznost aproksimacije. V tem prispevku uporabimo naslednjo obliko ta utežne funkcije: Weighting functions wk(x) must, in any direction, monotonously decrease with the distance from the corresponding sampling points xk and must be have a certain degree of continuity, which ensures continuity of the approximation. In the present work, we use the same form for all the weighting functions: n k=\ k=\ Uporaba metode premičnih najmanjših kvadratov - The Use of Moving Least Squares 585 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 (x,d) Nf(x-xk)i d, (12). Parametri d. v enačbi (12) določajo dejanski doseg vpliva vzorcev v različnih koordinatnih smereh. Različne oblike w(r) lahko pripravno prilagodimo namenu. Primera sta Gaussova in nasprotna polinomska oblika (si. 1): The parameters di in Equation (11) define the effective influence range of samples in the coordinate directions. A variety of forms of w(r) can be conveniently adapted, depending on the purpose. Examples are the Gaussian and reciprocal polynomial forms (Figure 1): wp(r) 1 1 + 1 r p , p = 2,3, 4,... (13). Mogoča je tudi uporaba funkcij s končno podporo, s čimer popolnoma izločimo vpliv oddaljenih točk na vrednost aproksimacije v dani točki. Temu se v glavnem izogibamo in za predstavljena področja uporabe uporabljamo funkcije, ki se dovolj hitro zmanjšujejo z razdaljo (na primer wG in wj, kar je po naših izkušnjah zadovoljivo. Takšne funkcije se dobro obnesejo pri neenakomerni porazdelitvi vzorčnih točk, ko na območjih z manjšo gostoto bolj oddaljene točke zagotovijo, da ostane sistem za določitev koeficientov (11) dobro definiran. Ker prostorska odvisnost uteži zagotavlja krajevno prilagajanje aproksimacije vzorčenim podatkom, ni potrebna velika množica osnovnih funkcij. Pri veliki gostoti vzorčnih točk in majhnih napakah lahko zadoščajo že monomi do prvega reda. Uporaba kvadratične polinomske osnove je Functions with a finite support can be utilized as well, which completely eliminates the distant points from affecting the approximation at the point of evaluation. We mainly avoid this and find functions with a sufficiently strong decay (e.g., wG and w4) adequate for the presented applications. These functions are handy for situations with a nonuniform distribution of samples, where in the region with small concentrations more distant samples ensure that the system for a determination of the coefficients (10) remains well posed. Since a spatial variation of weights ensures a local accommodation of the approximation to the sampled data, a large set of basis functions is not necessary. When the sample density is large and the errors are small, a linear basis consisting of monomials up to the first order may be sufficient. Experience shows that the quadratic polynomial basis, which we SI. 1. Obliki utežnih funkcij wjr) in w/r) Fig. 1. Weighting function forms wjr) and w/r) 586 Grešovnik I. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 glede na izkušnje primerna za veliko število različnih problemov. Dejanski doseg vpliva je odločilen parameter aproksimacije, katerega izbira je tesno povezana z izbiro osnove in lastnostmi aproksimirane funkcije. Enostavno vodilo pri tem je, da doseg vpliva ne sme biti večji od velikosti območja, na katerem lahko linearna kombinacija osnovnih funkcij z nespremenljivimi koeficienti dobro aproksimira funkcijo. Po drugi strani mora biti dejanski doseg v prisotnosti šuma čim večji, da se pri aproksimaciji izravna vpliv naključnih napak. SI. 2 shematično prikazuje metodo premičnih najmanjših kvadratov, pri katerih apsoksimiramo osem vrednosti s polinomskimi baznimi funkcijami drugega reda {1, x in x2}. Utežne funkcije pripadajoče vzorčnim točkam so narisane v spodnjem delu slike. Za izbrano točko izračuna aproksimacije so označene vrednosti uteži za vplivne vzorčne točke. used in the presented examples, is suitable for a large range of problems. The effective influence range is a crucial parameter for the approximation, and it is inseparably related to the chosen basis and properties of the approximated function. In simple terms, the range of influence should not be larger than the range on which a linear combination of basis functions with constant coefficients can adequately approximate the function. In the presence of noise, on the other hand, the effective range must be as large as possible in order to level out the effect of random errors. Figure 2 schematically presents the moving least-squares method where a set of eight sampled points are approximated with quadratic polynomial basis functions {1, x and x2}. The weighting functions corresponding to the sampling points are plotted in the lower part of the figure. For a chosen evaluation point, the values of the weights corresponding to the influencing samples are indicated. 2 ŠTUDIJA: APROKSIMACIJA ANALITIČNE FUNKCIJE 2 STUDY: APPROXIMATION OF AN ANALYTICAL FUNCTION V tem poglavju je obravnavana uporaba aproksimacije po postopku premičnih najmanjših In this section the application of the moving least-squares approximation is demonstrated on SI. 2. Shematičen prikaz aproksimacije premičnih najmanjših kvadratov Fig. 2. Scheme of the moving least-squares approximation Uporaba metode premičnih najmanjših kvadratov - The Use of Moving Least Squares 587 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 kvadratov na podlagi podatkov, ki so dobljeni z vzorčenjem analitične funkcije z dodanim naključno napako: /(jt)=sin(4x) + 0,5e0'- kjer sta rnd() enakomerno porazdeljena naključna napaka na območju [0,1] in R stalnica, ki določa raven šuma. Funkcijo vzorčimo na območju [Xj, xr]=[0, 5] v različnih številih m enakomerno oddaljenih točk. Proučujemo lastnosti aproksimacije glede na raven šuma R, število vzorčnih točk m in dejanski vplivni doseg d. Uporabimo kvadratično polinomsko osnovo: data sampled from an analytical function with an added random error: 2R(rndQ-0,5 (14), where rnd() is a uniformly distributed random number on [0,1] and R is a constant that defines the level of noise. The function is sampled over the interval [xl, xr]=[0, 5] using different numbers of equidistant points m. The properties of the moving least-squares approximation are studied with respect to the level of noise R, the number of sampling points m and the effective influence range d. An approximation with a quadratic polynomial basis is studied, i.e.: f1(x) = l,f2(x) = X,f3(x) (15). Za primerjavo rezultatov uvedemo mero za normalizirano gostoto vzorčenja, ki določa število vzorčnih točk na dejanski doseg d: p = m Kot mero za aproksimacijo napake uporabimo koren iz povprečne vrednosti vsote kvadratov odstopanj aproksimacije od/v Af =600 točkah: For a comparison of the results, a measure of the normalized sampling density is used that specifies the number of sampling points per effective influence range d: (16). d As a measure of the approximation error, the root mean square of the deviation between f and its approximation, calculated over a set of Ne=600 equidistant evaluation points, is used: l(yM-fM)1 X. = x, +(i — l)!L 1 J N 1 (17). Najprej študiramo aproksimacijo brez prisotnosti šuma. Aproksimacijo izračunamo za različna števila vzorčnih točk, pri katerih izberemo dejanski doseg tako, da je normirana gostota vzorčenja stalna, in sicer p — 1. SI. 3 prikazuje aproksimacijo z 10 in 15 točkami, ki ju primerjamo z aproksimirano funkcijo. Aproksimacijske napake so za različna števila točk prikazana v preglednici, ki je priložena k SI. 3, ter na SI. 4. Ko presežemo določeno število točk, se napaka enakomerno zmanjšuje z m. Tako je zaradi tega, ker se dejanski doseg zmanjšuje z naraščajočim m (saj se p ne spreminja) in ker lahko s končno polinomsko osnovo bolje aproksimiramo funkcijo na manjšem območju. V nadaljevanju obravnavamo aproksimacijo na podlagi podatkov s šumom. SI. 5 prikazuje aproksimacijo na podlagi 20 vrednosti z ravnijo šuma R=0A pri treh različnih vrednostih dejanskega dosega d. Jasno je viden vpliv dejanskega dosega An approximation without the presence of noise is studied first. This approximation is calculated for different numbers of samples where the effective influence range is chosen in such a way that the normalized sampling density is kept constant at p = 1. Figure 3 shows approximations with 10 and 15 points compared to the approximated function. Approximation errors using different numbers of sampling points are shown in the table that is attached to Figure 3 and in Figure 4. After some number of sampling pints m, the error monotonously falls with m. This is because the effective influence range falls with increasing m (since p is kept constant) and the finite polynomial approximation basis can be a better approximated function over a smaller interval. In what follows an approximation based of noisy sampled data is studied. Figure 5 shows approximations based on 20 sampled values with a noise level of R=0A, for three different effective influence ranges d. The impact of d on the 588 Grešovnik I. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 SI. 3. Aproksimacija na podlagi vzorčenih vrednosti brez šuma. Izmerki so narisani le za primer z 10 točkami, aproksimacija pa za 10 in 15 točk. Fig. 3. Approximation of sampled values without noise. Measurements are plotted only for the case with 10 sampled values, while the approximation is plotted for 10 and 15 sampled values. SI. 4. Odvisnost napake aproksimacije O od števila vzorčnih točk za vzorčenje brez šuma in pri stalnem p =1,0. Dodan je bolj podroben prikaz pri večjih m. Fig. 4. Dependence of the approximation error O on the number of sampling points for sampling without noise, at constant p =1.0, with a detailed plot for larger m. na kakovost aproksimacije. Ko je d prevelik, postanejo izbrane osnovne funkcije nezadostne za aproksimacijo funkcije na območjih velikostnega reda 2d okrog točke izračuna, na katerih vzorčne točke pomembno prispevajo k aproksimaciji (enačbe (11), (12), (13)) Na SI. 5 se to kaže v preveč zravnani aproksimaciji pri največjem d=0,4, ki ne zmore slediti nihanju aproksimirane funkcije. Ko je d premajhen, se aproksimacija nagiba k approximation quality is clearly visible. When d is too large, the approximation basis becomes insufficient for adequately approximating the function over intervals of order of magnitude of 2d, on which the sample contributions are significant (consider Equations (10), (11), (12)). In Figure 5 this is reflected in a somehow flattened approximation with the largest d=0.4, which cannot follow the oscillations of the approximated Uporaba metode premičnih najmanjših kvadratov - The Use of Moving Least Squares 589 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 Sl. 5. Učinek dejanskega dosega vpliva vzorčnih točk na aproksimacijo na osnovi 20 vzorčnih točk s šumom Fig. 5. Impact of the effective influence range of samples on an approximation based on 20 sampled values with noise interpolaciji vzorčenih vrednosti in zato sledi tudi naključnim spremembam, ki so posledice šuma. Videti je, da pri dani funkciji, ki jo aproksimiramo, izbranih vzorčnih točkah in določeni ravni šuma obstaja optimalen dejanski doseg, pri katerem je kakovost aproksimacije največja. To je vidno tudi na sliki 7, kjer je prikazana odvisnost mere za napako aproksimacije a od function. In contrast, when d is too small, the approximation tends to interpolate the sampled data and therefore follows the spurious fluctuations caused by the random noise. It seems obvious that for a given approximated function, sampling points and level of noise, there exists an optimal effective range for which the approximation quality is the best. This is Sl. 6. Aproksimacija na osnovi 40 vzorčnih točk z višjo ravnijo šuma (R=0,8) Fig. 6. Approximation based on 40 samples with a larger noise level (R=0.8) 590 Grešovnik I. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 povečanim detajlom na desni strani slike Fig. 7. Dependence of the error measure O on the reciprocal effective influence range for data from Figure 6, with enlarged detail on the right-hand side količine, ki je obratno sorazmerna z d. Graf se seen from Figure 7, where the dependence of the nanaša na višjo raven šuma (L=0,8) in število error measure a on reciprocal d is plotted. The vzorčnih točk (m=40). Za vsako točko grafa so graph refers to a larger level of noise (L=0.8) and uporabljene iste vzorčene vrednosti. Praviloma bi number of sampling points (m=40). The same morali vrednosti na grafu povprečiti po več sampled data is used for each point in the graph. naključnih izidih vzorčenih vrednostih. Graf bi v Otherwise, the values on the graph should be tem primeru vseboval naključne spremembe, averaged over a number of sampling realizations vendar bi bila usmeritev podobna. SI. 6 prikazuje and the graph would contain random fluctuations, vzorčene vrednosti uporabljene pri sliki 7 in but the trend would be similar. Figure 6 shows the aproksimacijo na podlagi teh podatkov z dejanskim sampled data used in Figure 7 and an approximation dosegom d, ki je blizu optimalnega. of this data with d that is close to the optimal value. Podobno kakor v primeru brez šuma Similar to the data without noise, it is pričakujemo, da lahko kakovost aproksimacije expected that the quality of the approximation can izboljšujemo z večanjem števila vzorčnih točk m. be improved by increasing the number of sampling Vendar šum preprečuje neomejeno hitro večanje points m. However, the noise prevents an unlimited natančnosti z večanjem m. To se zgodi v območju, rapid improvement of the approximation accuracy kjer amplituda šuma doseže podoben velikostni red with increasing m. This happens at the point where kot povprečna aproksimacijska napaka v primeru, the amplitude of the noise reaches a similar ko ni šuma. Od tu naprej lahko aproksimacijsko magnitude to the approximation error for the case napako še vedno manjšamo z večanjem gostote without noise. From this point on the approximation vzorčenja, vendar je to posledica dejstva, da error can still be reduced by increasing the sampling naključne napake zaradi šuma povprečimo po density, but this is attributed to the fact that the večjem številu vzorčnih točk, ki jih zajamemo random sampling errors are averaged over a larger znotraj območja vpliva okrog vsake točke, v kateri number of sampled values, captured within the area izračunamo aproksimacijo. of influence around each calculation point. Opisan učinek je jasno viden na sliki 8, kjer The effects described above are clearly je narisana napaka aproksimacije pri različnih visible in Figure 8, where the approximation error številih vzorčnih točk. Do približno m=l 5 se napaka is plotted for different numbers of sampling points. hitro zmanjšuje z rastočim m zarasi zmanjševanja Up to about m=15 the error is rapidly decreasing razdalje med točkami. Učinek šuma v tem delu ni with growing m because of the decreased spacing izrazit, ker je premajhna količina podatkov between the sampling points. The effect of noise is poglaviten vir nenatančnosti. Pri večjih m je not clearly visible in this regime because insufficient Uporaba metode premičnih najmanjših kvadratov - The Use of Moving Least Squares 591 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 SI. 8. Aproksimacijska napaka v odvisnosti od števila vzorčnih točk pri R=0,2 in z nespremenljivim dejanskim dosegom vpliva d=5/17. Vsaka točka na grafu se nanaša na en sam izid zbiranja vrednosti. Fig. 8. Approximation error dependent on the number of measurements for data for R=0.2 and with a constant effective influence range d=5/17. Each point on the graphs refers to a single realization of the sampled data. nadaljnje izboljševanje natančnosti omejeno z naključnimi napakami v podatkih. Aproksimacijske napake se počasi zmanjšujejo s povečevanjem gostote vzorčenja in naključno nihajo zaradi stohastične narave vzorčenja. Na SI. 8 vsaka točka na grafu predstavlja en sam izid zajemanja podatkov, ki vsebujejo naključne napake. Usmeritev počasnega povprečenja napak pri povečani gostoti vzorčenja je vseeno razvidna, ker je prikazano večje število točk. data is the main cause of the inaccuracy. At larger m, a further increase of accuracy is restricted by the random errors in the sampled data. Approximation errors tend to slowly decrease with the increasing sampling density and fluctuate randomly because of the stochastic nature of the sampling. In Figure 8 each point on the graph corresponds to a single realization of sampled data containing random noise, and the average trend for the error can be observed because of the large number of presented points. 3 PRIMERA UPORABE 3 EXAMPLE APPLICATIONS 3.1 Glajenje meritev pri poskusu 3.1 Smoothing of the Experimental Measurements Na sliki 9 je prikazana uporaba aproksimacije po metodi premičnih najmanjših kvadratov pri glajenju signala zaznavala za merjenje temperature. Podatki so iz [5] in so bili uporabljeni za obrnjeno določitev parametrov prestopa toplote in trenja med kovinami z mazanjem iz meritev pri laboratorijskem testu zoževanja traku. Test je bil simuliran z metodo končnih elementov, iskani modelski parametri pa so bili določeni z minimizacijo mere za odstopanje meritev od simuliranih podatkov, definirane kot: Figure 9 shows the application of the moving least-squares approximation for smoothing the output from a temperature sensor. The data is from [5] and was used for an inverse reconstruction of the heat-transfer and friction parameters between lubricated metals from the measurements performed with the strip-reduction laboratory test. The laboratory test was simulated using the finite element method and the unknown parameters were obtained by minimizing a measure of the discrepancy between measured and simulated data, defined as: f ( p ) =XJttmax Wi ( Mim t-Mi c p,t )) dt (18). 592 Grešovnik I. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 V zgornji enačbi i teče po uporabljenih časovno odvisnih merjenih količinah, Mi m (t) je ustrezna izmerjena količina (sila na trak ali temperatura v dani točki), M^ je ustrezna količina, izračunana s simulacijo preizkusa pri preizkusnih vrednostih iskanih parametrov, in W. uteži, izbrane za posamezne vrste meritev. Ker je bila pri numerični simulaciji uporabljena integracija s prilagodljivim korakom in zaradi merskega šuma (si. 9) vsebuje funkcija Xp) naključna nihanja, kar preprečuje učinkovito uporabo numeričnih minimizacijskih metod. Ta problem je bil odpravljen z zamenjavo interpoliranih meritev z gladkimi aproksimacijami po metodi premičnih najmanjših kvadratov (nepretrgana črta na SI. 9). Pri aproksimaciji je bila uporabljena kvadratična onova (15) z utežnimi funkcijami oblike wjr) iz (13). Dejanski doseg vpliva je bil postavljen na d = 0,1 s. Na ta način nismo zameglili prehodnih pojavov pri časovnem poteku meritev, hkrati pa je bil vpliv šuma dovolj izravnan (povečan detajl na sliki 9), da smo dobili gladko odzivno funkcijo /, pri kateri smo lahko učinkovito določili minimum. In the above equation, i runs over the measured quantities, Mi(m)(t) is the corresponding measured quantity (force on the strip or the temperature at a certain point), Mi(c) is the corresponding quantity obtained by the finite-element simulation with the trial values of the searched parameters p and Wi weights assigned to different kinds of measured data. Because an adaptive time-integration scheme was applied in the numerical simulation and due to the measurement noise (Figure 9), the minimized function f(p) contained random oscillations that prevented an efficient numerical minimization. This problem was alleviated by replacing the interpolated measurement data with a smooth moving least-squares approximation (the solid line in Figure 9). A quadratic basis (14) was used with the weighting function wG(r) from (12). The effective influence range was set to d = 0,1. In this way the transient features of the sensor response were not smeared, while the effect of the noise was sufficiently leveled (the enlarged detail in Figure 9) to obtain a smoothed response function f that could be efficiently minimized. 3.2 Aproksimacija odzivnih funkcij pri optimizaciji 3.2 Approximation of the Response Functions During Optimization V naslednjem primeru obravnavamo optimizacijo, pri kateri imamo opravka z odzivnimi funkcijami s šumom. Podatki (SI. 10) so iz [6], kjer smo optimirali širino in višino kanala, ki je izdelan The next example treats optimization with noisy response functions. The data (Figure 10) is taken from [6], where the width d and the height h of a channel produced by blow forming was SI. 9. Zglajen signal (nepretrgana črta) temperaturnega zaznavala s povečanim detajlom na desni strani Fig. 9. Smoothed data (solid line) from the temperature sensor with enlarged detail on the right-hand side Uporaba metode premičnih najmanjših kvadratov - The Use of Moving Least Squares 593 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 a) b) SI. 10. Namenska (a) in omejitvena funkcija (b) optimizacijskega problema, izračunani z metodo končnih elementov na pravilni mreži 20x20 točk Fig. 10. Objective (a) and constraint function (b) of the optimization problem calculated using a FEM simulation on a regular 20x20grid z napihovanjem, z namenom zmanjšati tveganje za nastanek razpok, ki se pojavijo zaradi lokalizacije deformacije. Naloga je bila postavljena kot maksimizacija površine prereza kanala S.(d,h) (SI. 10 a)) pri omejitvi, da je tlak P.(d,h) , pri katerem se začne lokalizacija, pod predpisano mejo. Zaradi tehničnih zahtev sta parametra omejena z 8 mm < d < 12 mm in 3 mm < h < 3,4 mm. Šum v odzivnih funkcijah izvira iz numerične simulacije postopka napihovanja, ki je uporabljena za izračun odzivnih funkcij [6]. Uporabiti je bilo treba prilagodljivo izboljševanje mreže končnih elementov z veliko gostoto elementov v območju lokalizacije, zaradi česar bi bilo težko zmanjšati raven šuma. Da smo lahko poiskali optimalno rešitev, smo najprej izračunali odzivni funkciji v pravilni mreži točk. Dobljene podatke smo zgladili z aproksimacijo po metodi premičnih najmanjših kvadratov z obliko utežnih funkcij wG(r), dejanskima dosegoma dl = 1 mm in d2 = 0,1 mm v smereh d in h ter kvadratnimi polinomskimi osnovnimi funkcijami: optimized in order to reduce the risk of localization-induced defects. The task was formulated as a maximization of the channel cross-section Si(d,h) (Figure 10 a)) using a constraint that the forming pressure Pi(d,h) at which localization occurs (Figure 10 b)) is below some prescribed limit. Due to technical requirements the parameters were limited to 8 mm ? d ? 12 mm and 3 mm ? h ? 3,4 mm. The noise originating from the finite-element numerical simulation of the blow-forming process was applied for the calculation of the response functions [6]. Adaptive mesh refinement with a high mesh density at the localized zone had to be applied, because of which it would be difficult to reduce the level of noise. In order to find the optimal solution, the response functions were first sampled on a regular grid of points. The sampled data was smoothed by the moving least-squares approximation with the weighting function wG(r), effective ranges d1 = 1 mm and d2 = 0,1 mm (corresponding to d and h) and the basis functions: f(x) = 1; f2(x) = xi; f3(x) = x2; f4(x) = xi2; f5(x) = x22; f(x) = xx (19), kjer sta x=d in x=h. Aproksimiran odziv je prikazan na slikill. V [6] je bil optimizacijski problem rešen z aproksimiranim odzivom z metodo zaporednega kvadratičnega programiranja. Tak postopek je where x1=d and x2=h. The approximated response is shown in Figure 11. In [6] the optimization problem with smooth approximated response functions was solved by a sequential programming algorithm. Such an 594 Grešovnik I. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 c) d) Sl. 11. Zglajene odzivne funkcije s slike 10 Fig. 11. Smoothed response functions from Figure 10 neučinkovit, ker zahteva dovolj gosto vzorčenje approach is inefficient because it requires a sufficiently odziva po celotnem dovoljenem območju v prostoru dense sampling of the response over the whole allowed parametrov. Pri večjem številu parametrov to range of parameters. If the number of parameters is postane zdaleč prezahtevno. large, this becomes prohibitively expensive. Zato smo uresničili drugačen iterativen A different iterative solution approach was, postopek, pri katerem funkcije aproksimiramo therefore, designed, where the response approximation krajevno na omejenem območju. V vsaki iteraciji is maintained locally in a limited region of interest. določimo omejeno območje zanimanja. Na For each iteration a restricted region of interest is podlagi vzorčenih vrednosti odzivnih funkcij, ki defined. The response functions are calculated smo jih izračunali v trenutni in predhodnih (sampled) at a number of points within the region of iteracijah, aproksimiramo odzivne funkcije po interest. Based on the sampled values from the current metodi premičnih najmanjših kvadratov. Ustrezno and previous iterations, the moving least-squares kakovost aproksimacije želimo doseči le na approximations of the response functions are trenutnem območju zanimanja, kar dosežemo z constructed. A good approximation quality is ustrezno nastavitvijo dejanskega dosega. Nato considered only within the current region of interest, rešimo optimizacijski problem, pri katerem which is achieved by appropriately setting the effective nadomestimo prvotno namensko in omejitvene range of influence. The optimization problem is then funkcije z ustreznimi aproksimacijami in dodamo solved, with the original objective and constraint omejitev koraka, s čimer omejimo možne rešitve functions being replaced by the approximated ones. na trenutno območje zanimanja. Rešitev tega The step restriction is also added, which restricts the notranjega optimizacijskega problema postane nov possible solutions to the current region of interest. The približek za rešitev prvotnega problema in središče solution of this internal problem becomes the new področja zanimanja v naslednji iteraciji. Velikost current guess and the center of the region of interest območja povečamo, če leži rešitev na robu in the next generation. The size of this region is območja in zmanjšamo, če leži dovolj vstran od increased if the solution lies on the edge of the region meje. Postopek ponavljamo, dokler ne dosežemo of interest, and decreased if it lies far enough from the predpisane natančnosti, ali dokler ni več mogoče edge. The procedure is repeated until the prescribed nadaljnje izboljšanje natančnosti glede na raven accuracy is reached or no further improvement is šuma. possible, with respect to the level of numerical noise. Konvergenca opisanega algoritma je The convergence of the described algorithm prikazana na sliki 12. Območja zanimanja in is shown in Figure 12. Regions of interest and the vzorčene vrednosti so prikazane po iteracijah. V sampled values through the iterations are indicated. Uporaba metode premičnih najmanjših kvadratov - The Use of Moving Least Squares 595 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 prvi iteraciji določi območje zanimanja uporabnik, For the first iteration the region of interest is defined število vzorčnih točk pa je izbrano tako, da je by the user and the number of samples is chosen in nekoliko večje od števila baznih funkcij. V such a way that it is slightly larger than the number of naslednjih iteracijah je število novih točk na basis functions. In the subsequent iterations, the iteracijo za eno večje od števila parametrov, v našem number of new sampling points per iteration is set to primeru 3. be one more than the number of parameters, i.e., 3. Vzorčne točke so naključno izbrane znotraj The sampling points are chosen randomly within območja zanimanja. To je mogoče, ker pri metodi the region of interest. This is possible because no premičnih najmanjših kvadratov ni potrebna particular arrangement of the sampling points is required kakšna posebna razporeditev točk. Treba je for the moving least-squares approximation. It must be omeniti, da so možne numerične težave povezane mentioned that numerical difficulties related to a badly s slabo pogojenim sistemom enačb (11) za izračun conditioned system of equations for the calculation of koeficientov aproksimacije. Sistem lahko postane approximation coefficients (10) could potentially occur. slabo pogojen, če je število vzorčnih točk z The system can become badly conditioned if either the znatnimi utežmi premajhno, ali če so točke number of sampling points with significant weights is razporejene na poseben način, tako da so osnovne insufficient or the points are arranged in a special way funkcije na množici točk skoraj linearno odvisne. such that the basis functions are close to being linearly Prvo možnost lahko učinkovito preprečimo z dependent on the set of points. The first situation is uporabo utežnih funkcij, ki se ne zmanjšujejo efficiently prevented by using weighting functions that prehitro z naraščajočo razdaljo od vzorčnih točk. do not decay too rapidly with the distance from the Tako ima vedno dovolj točk iz prejšnjih iteracij sampling points. In this way enough sampling points povsod tam, kjer računamo aproksimacijo, dovolj from previous iterations will always have large enough SI. 12. Konvergenca optimizacijskega algoritma, ki temelji na zaporednih aproksimacijah odzivnih funkcij. Pravokotniki označujejo območja vzorčenja v zaporednih iteracijah, krožci prikazujejo točke, v katerih so bile izračunane odzivne funkcije, večji krožci pa rešitve aproksimiranih problemov. Za ponazoritev so prikazane tudi obrisi zglajenih odzivnih funkcij (SI. 10). Fig. 12. Convergence of the optimization algorithm based on successive response approximations. Rectangles denote sampling regions within successive iterations, dots represent points where the response was sampled and larger dots represent the minima of successive approximated problems. The contours of smoothed response (Figure 10) are plotted for orientation. 596 Grešovnik I. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 velike vrednosti utežnih funkcij, da se izognemo nezadostnosti podatkov. Glede na izkušnje je w4(r) iz (13) primerna oblika utežnih funkcij. Druga možnost je zelo malo verjetna, kadar uporabljamo naključno vzorčenje, posebej, če je število vzorčnih točk z znatnimi utežmi precej večje od števila osnovnih funkcij. V predstavljenem primeru smo dobili rešitev v 10 iteracijah, v katerih sta bili odzivni funkciji izračunani v 36 točkah. To je precej bolj učinkovito kakor izvedba optimizacije na aproksimiranem odzivu po celotnem dovoljenem območju. Pri večjem številu parametrov postane razlika v učinkovitosti izrazita in aproksimacija po celotnem območju neizvedljiva. 4 SKLEP V prispevku je nakazano, da je metoda premičnih najmanjših kvadratov vsestranska in prilagodljiva aproksimacijska metoda, ki je zaradi posebnih značilnosti uporabna pri reševanju številnih praktičnih problemov. Za aproksimacijo ni potrebna kakšna posebna ureditev vzorčnih točk ali delitev območja aproksimacije. Razmeroma majhno število osnovnih funkcij zadošča za dovolj natančno aproksimacijo gladkih funkcij na poljubno velikem območju. Velikost sistema enačb za določitev koeficientov aproksimacije se ne poveča, če povečamo gostoto vzorčenja, zato pa moramo sistem enačb rešiti posebej v vsaki točki, v kateri izračunamo aproksimacijo. Povečevanje natančnosti z gostejšim vzorčenjem ne prizadene stabilnosti aproksimacije. Metoda je zato sorazmerno bolj učinkovita, kadar je potrebna večja natančnost, aproksimacijo pa je treba izračunati v manjšem številu točk. Zveznost aproksimacije je odvisna od zveznosti osnovnih in utežnih funkcij. Kadar so oboje gladke, so tudi višji odvodi zvezni. Z uporabo utežnih funkcij, ki se počasneje zmanjšujejo z razdaljo, lahko dosežemo stabilnost metode tudi v primerih, pri katerih ne moremo zagotoviti enakomerne gostote vzorčenja po celotnem območju aproksimacije. Pri določeni izbiri baznih funkcij se filtrirajo višjefrekvenčne oscilacije glede na dejanski doseg vpliva vzorcev, kar lahko uporabimo za izravnavo vpliva šuma v vzorčenih podatkih. Večji dejanski doseg pomeni boljše glajenje, vendar tudi slabšo zmožnost prilagajanja prehodnim lastnostim weights at all points of the evaluation in order to prevent data deficiency. Based on experience, w4(r) from (12) is a suitable form for the weighting function. The second situation is very unlikely when random sampling is used, especially when the number of sampling points with significant influence is considerably larger than the number of basis functions. In the presented case the algorithm converged in 10 iterations in which 36 evaluations of the response were performed. This is significantly more efficient than performing an optimization of the approximated response over the whole permitted domain. With a large number of parameters, the difference in the efficiency becomes drastic and the approximation over larger domains becomes unfeasible. 4 CONCLUSION We have shown that the moving least-squares method is a versatile approximation technique suitable for many practical applications because of its distinctive set of features. No regular grid of points or partition of the domain is necessary. A relatively small number of basis functions can be used for an accurate approximation of smooth functions over arbitrarily large domains. The size of the system of equations for the approximation coefficients is not increased when increasing the sampling density, at the cost that the system must be solved for every evaluation point. The stability of the approximation is not affected when improving the accuracy by increasing the sampling density. The method is therefore relatively more efficient when higher accuracy is required, but the approximation is evaluated fewer times. The continuity of the approximation depends on the continuity of the basis and weighting functions. When both are smooth, the approximation does not have discontinuous higher derivatives. By using weighting functions with a slower decay, good stability of the method is achieved in cases when a uniform sampling density cannot be ensured over the approximation domain. For a given set of basis functions, higher-frequency oscillations are filtered, depending on the effective influence range of the samples, which can be used to compensate for the effects of noise in the data. A larger effective range means better smoothing, but also a poorer ability to follow the transient features of the approximated function. By Uporaba metode premičnih najmanjših kvadratov - The Use of Moving Least Squares 597 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)9, 582-598 aproksimirane funkcije. S postavitvijo dejanskega dosega na red velikosti razdalje med vzorci ali manj lahko dosežemo približno interpolacijo podatkov. V praksi moramo doseči primeren kompromis med opisanimi učinki z ustrezno nastavitvijo dejanskega dosega glede na raven šuma, gostoto vzorčenja in lastnosti aproksimirane funkcije. Aproksimacijo po metodi premičnih najmanjših kvadratov lahko preprosto izboljšujemo z dodajanjem vzorčnih točk. Zaradi tega je metoda posebej primerna za uporabo v optimizacijskih postopkih, ki temeljijo na zaporedni aproksimaciji odzivnih funkcij setting the effective influence range at the order of the spacing between the samples and below, the approximation tends to interpolate the data. In practice, a suitable compromise between these effects must be achieved by appropriately adjusting the effective influence range with respect to the level of noise, sampling density and properties of the approximated function. The moving least-squares approximation can be easily refined by the addition of new sampling points. This makes it particularly suitable for use in optimization methods with a successive adaptive approximation of the response functions. 5 LITERATURA 5 REFERENCES [1] W.H. Press, B.P. Flannery, SA. Teukolsky, W.T. Vetterling (1992) Numerical recipes in C: The art of scientific computing. Second Edition, Cambridge University Press. [2] I. N. Bronstein, K. A. Semendjajew, G. Musiol and H. Muehlig (2000) Handbook of mathematics, 5th Edition. Harri Deutsch. [3] P. Breitkopf, A. Rassineux, P. Villon (2002) An introduction to the moving least squares meshfree methods, Revue Europenne des lments Finis, Volume 11 - No 7-8, 2002. [4] I. Emri, R. Cvelbar (2006) Uporaba gladilnih funkcij za glajenje podatkov podanih v diskretni obliki. Strojniški vestnik - Journal of Mechanical Engineering, Vol. 52, p.p. 181-194, 2006. [5] I. Grešovnik, S. Stupkiewicz, T Rodič (2006) Sensitivity analysis and inverse modelling of limits of lubrication. FP6 European project ENLUB (G1RD-CT-2002-00740), Development of New Environmentally Acceptable Lubricants, Tribological Tests and Models for European Sheet Forming Industry. Final technical report. [6] I. Grešovnik, S. Hartman, T. Rodič (2005) Reducing localisation induced defects at blow forming in terms of optimal shape design. In E. Onate (ed.), Computational plasticity: fundamentals and applications, Proceedings of the eighth international conference on computational plasticity held in Barcelona, Spain, 5th-7th September, 2005. Part 1, p.p. 388-391. [7] R. Fletcher (1996) Practical methods of optimization (second edition). John Wiley & Sons, New York. Avtorjev naslov: dr. Igor Grešovnik Center za računalništvo v mehaniki kontinuuma - C3M p.p. 431 1102 Ljubljana igor@c3m.si Author’s Address: Dr. Igor Grešovnik Centre for Computational Continuum Mechanics - C3M P.O. Box 431 1102 Ljubljana, Slovenia igor@c3m.si Prejeto: 16.9.2006 Received: Sprejeto: 25.4.2007 Accepted: Odprto za diskusijo: 1 leto Open for discussion: 1 year 598 Grešovnik I.