Új módszer kísérletek optimális tervezésére

(1)

ÚJ MÓDSZER

KlSERLETEK OPTIMÁLIS TERVEZÉSÉRE

DR. SZIDAROVSZKY FERENC

A gyakorlati élet számos területén van szükség arra, hogy véges sok pontban ismert függvény újabb pontbeli értékét vagy egy adott tartományon való átlagértékét megbecsül—

jük. Például az ásványvagyon mennyiségi és minőségi paramétereinek becslése is ezt afelada- tot jelenti, ahol a véges sok pontbeli függvényértéket a kutatófúrások adataival azonosítjuk.

Minthogy újabb függvényérték meghatározása általában igen költséges, fontos feladatot jelent a kutatás olyan jellegű irányítása, hogy minél kevesebb függvényérték felhasználásá- val a lehető legtöbb információt nyerjük. Ezt az utóbbi feladatot úgy is megfogalmazhatjuk, hogy a becsült eredmény bizonytalanságát minél kevesebb függvényérték felhasználásával a lehető legkisebbre szorítsuk. Természetesen felvetődik e kérdésnek a következő meg—

fordítása: hogyan, milyen kisérlet, kutatástervezés mellett juthatunk minimális, de alkal—

masan választott mérési hely felhasználásával előirt bizonytalansági szint alá. '

Ebben a dolgozatban az itt felvetett kérdésekre adunk választ olyan feltételek mellett, hogy becslési módszerként a statisztikai irodalomból ismert Kriging becslési eljárást vá- lasztjuk. Mielőtt a részletek ismertetésére sor kerül, a Kriging módszer alapjait vázoljuk.

'1. Kriging módszer alapelve

jelöljük Z-vel azt a függvényt, amelynek átlagértékét kívánjuk meghatározni egy T tartományon, ahol T az m dimenziós tér korlátos, zárt halmaza, és Z integrálható T-n.x Tegyük fel, hogy Z értékét a diszkrét x1, xz, . . ., x,, pontokban megmértük, és jelölje Zl, ZZ, ., ., Z,, a mért függvényértékeket. Legyen x valamilyen hely Z értelmezési tarto- mányában. A Z(x) függvényértéket úgy tekintjük, mint egy sztochasztikus függvény x

helyen való realizálását, amelyre feltesszük (lásd (1), (2), (3)), hogy

E

(

Z x

( ))

: M,

/1/

Cov (Z(x-l—h), Z(x)) : K(h),

ahol E a várható értéket, Cov a kovarianciát, h pedig h hosszát jelöli. Az /1/ összefüggésben

tehát azt tesszükfel, hogy Z(x) várható értéke független x-től, valamint a Z(x—l—h), Z(x) értékek kovarianciája csak az abszcisszák távolságától függ. Legyen továbbá

win : % Var (Z(x—I—h)—Z(X)), /2/

ahol Var a variancia jele. Könnyen belátható, hogy az (1) feltételek mellett

y(h) : K(0)——K(h), /3/

(2)

DR. SZIDAROVSZKY: KISÉRLETEK TERVEZÉSE 1115

valamint y(h) értéke konkrét h : 0 mellett a

N(h)

van 2501)_É (learn—zal)? /4/

formulával közelíthető, ahol N(h) jelöli azon mérési pontpárok számát, amelyek távolsága éppen h. A /4/ becsléssel a konkrét alappontokban meghatározotty függvényt a legkisebb négyzetek módszerével ezután simítjuk. A gyakorlatban legtöbbször alkalmazott, így ott ún. varíogramtípusok:

y(h) : wha (0 : oz ( 2) ' (polinomiálls variogram) /5/

y(h) : con—eü) (exponenciális variogram) /6/

'y(h) : w(1—-e""*") (Gauss típusú variogram) /7/

3

j: (ÉL—(i)) ha h§oz

))(h) : 2 a az , (szférikus variogram) /8/

lm, ha h ) a

Megjegyezzük, hogy a geostatisztikában az ásványvagyon becslésekor általában a /8/ vario- gramtípust alkalmazzák.

Vezessük be ezután a következő jelöléseket:

vü- -—— Mllxi—lel), m— :;j—íj mux —xn)dx,

" 1 1

ym : Tf— J Jyle—x'll) dx' dx, Z: _lfl— fZ(x) dx, /9/

T T T

területét, térfo—

gatát egy-, két-, illetve háromdimenziós esetben) jelöli. Keressük a Z függvény T-n vett átlagértékének becslését:

Z : Z zizi /10/

ia).

alakban, ahol a I'll, . . ., Á,, konstansok meghatározása jelenti a becslési eljárást. Feltesszük, hogy a /10/ becslés torzítatlan, azaz

M : Z AM, /11/

i:1 .

amely nyilvánvalóan fennáll, ha

2 2,- : 1. /12/

isl

A /11/ összefüggésben felhasználtuk, hogy 2 várható értéke is M, azaz a várható érték

képzés és a T-n vett integrál felcserélhető. Ez pedig elég sima függvények esetén mindig teljesül. Feltételezzük, hogy a /10/ becslés minimális szórású is, azaz

E [(í—í; 24202] : min /13/

(3)

1 1 16 DR. SZlDAROVSZV FERENC

A /13/ célfüggvénynek a /12/ feltétel melletti minimalizálását a Lagrange-féle multipli—

kátor elv (4) felhasználásával visszavezethetjük a feltétel nélküli

E [li—Él Min—Zu [; li—1] _, min /14/'

szélsőérték-feladatra. Elvégezve az első tagbeli négyzetre emelést és a várható értéket tagonként képezve, valamint felhasználva a /3/ összefüggést, a célfüggvényre nyerjük:

n n ". n

_?700'l'2 ;: 273711?" 2 E ÁiÁj'ij—"ZN [Z Áj_1jl I'lS/x'

:1 izl 521 f:]-

amelynek a ,u és a 1,- változók szerinti parciális deriváltjait képezve a

/16f

"'n; §;st gli—m (5 : 1, 2. .. ., n)

lineáris egyenletrendszert kapjuk az ismeretlen Zj paraméterekre. Az egyenletrendszer"

tömör felírása érdekében vezessük be a

0 1 1 1 1 1 ,u,

1 0 712 ' ' ' 71,n—1 %n 3710 A1

1 721 0 ' - - Van—1 Van 720 Az

Gn : --- §?" : ... ' A" : . s . . /17/

1 yn—IJ lő.—14 - - - 0 lán—1." 77144) Árt—1

1 %u %a - - - ?)",n—l 0 %,o 2-7.

szimmetrikus matrixot és vektorokat. Ekkor a /16/ egyenletrendszer tömören

Gui." : vu /18/

alakban írható fel. Könnyen kimutathatjuk, hogy ha a)," vektort /18/ megoldásának választ-

juk, akkor a /13/ variancia-érték a lényegesen egyszerűbb

"* ; Win—?m : 343651 Ja.—7700 /19/

alakban is felírható, ahol most T—vel a transzponálást jelöltük (5). A gyakorlati alkalmazá- sokban a /18/ egyenletrendszert alkalmasan választott numerikus eljárás segítségével oldják

meg ((ő), (7)), és a megoldásként kapott A,— értékek felhasználásával a /10/ formula adja

Z átlagértékének becslését és /19/ a becslés varianciáját.

2. Adott számú mérési hely optimális tervezése

Helyezzünk a vizsgált tartományra elég sűrű hálót, és jelölje y1,y2, .. ., yN a csúcs- pontokat. Természetesen feltesszük, hogy a már meglevő mérési helyek e csúcspontok között szerepelnek. Ha előzetesen meglevő információink alapján valamilyen részhalmazra nem kívánunk újabb mérési helyet telepíteni, akkor az ide eső csomópontokat egyszerűen kihagyjukafelsorolásból. Legyen ezután GN, yN, AN a /17/ definicióknak megfelelő mennyi- ség, ha az X; alappontok helyett most az yj (] : 1, 2, .. ., N) pontokat szerepeltetjük.

(4)

KlSÉRLETEK TERVEZÉSE

1 l 17

Tegyük fel ezután, hogy r számú újabb mérési helyet kívánunk oly módon meghatá- frozni, hogy az ezek figyelembevételével nyert (19) típusú variancia-érték minimális legyen.

igen jó közelítését kapjuk a feladatnak, ha az újabb r pontot a még ki nem jelölt y,,jtl,

*ynjtg, . . ., yN csúcspontok közül választjuk. Az újonnan kiválasztott pontokat célszerűen egy 0 vagy 1 értékű komponensekből álló d : (do, dj, . . ., d", d,,H, . . ., dN) vektorral jellemezhetjük, ahol d,— : 1 (O § i § n), valamint egy továbbidi komponens 1 volta az i-edik pont kiválasztását jelenti, a dí : 0 eset pedig annak felel meg, hogy az i-edik pont nem kerül kiválasztásra. Ha pontosan r számú új pontot akarunk kijelölni, akkor az itt bevezetett Vd,- számokra nyilvánvalóan fennáll, hogy

d,:1 (i:1,...,n) d,e(0,1)(i:n-l-1,---,N) /20/

A7

í:;i1

d-:r._Z

Vegyük továbbá azt is észre, hogy ha egy d,- komponens 0, akkor a nekik megfelelő 'sort és oszlopot a GN matrixból, a megfelelő elemet a W és AN vektorból ki kell hagynunk, hiszen a Kriging módszer /16/ egyenletrendszerében csak a valóban bevont pontokhoz tar- tozó értékek szerepelhetnek. Minthogy az együttható matrix méretét nem célszerű meg- változtatnunk, sorok, oszlopok és elernek kihagyása helyett egyszerűen O—val való helyette- sítésüket célszerű elvégeznünk, hiszen ez csak triviálisan teljesülő 0 : O alakú egyenle- tekkel történő bővítést jelent. A megfelelő sorok, oszlopok és elemek zérussá való transz—

formációját a megfelelő di értékekkel való szorzással végezhetjük el leginkább. A d,— érté- kekkel való szorzás természetesen nemcsak a d,— : 0 esetben megfelelő, hanem a d,- : 1 esetben is, hiszen ekkor a di-vel való szorzás semmilyen változtatást sem jelent. Az elmon- dottak szerint tehát a /17/ egyenletnek az újabb pontok bevonása utáni alakja:

N

; ajd, :1

N F] /21/

dia: §; (y,, d,d,)(1,d,.) : %d,. (i : 1, 2, . . ., N),

,:

valamint a /19/ variancia is a hasonló

N

M- ;; (z,- d,-) (ne dió—700 /22/

alakban írható fel. Ily módon a d,- (] : n—H, .. ., N) paraméterek optimális megválasztása a következő matematikai programozási feladat megoldását igényli:

d,:1(í:1,....n) d,e(o,1) (i:n-j—1,...,N)

N

d ::

f:;a-l ] r

/23/

N

§: Aja, : 1

;:1

N

Wii; y,,d,d,i, : awm (i : 1, 2, . . ., N)

N

——')700—j—u-l— ,; A,. dia—% : min.

(5)

1118 DR. SZIDAROVSZY FERENC

Itt a d,- (i : 1, . . ., N), M és a Á,— (j : 1, . . ., N) számok megválasztása jelenti az isme—

retleneket. A Vi], f),-0 és 7700 mennyiséget ismertnek tételezhetjük fel, his2en azok a vario—-

gram-függvény alapján könnyen származtathatók. Feladatunk tehát egy vegyes programo—

zási probléma, ahol a d,- változók diszkrétek, 0 vagy 1 értékűek, a ,a és 2.) változók Pedig valósak. A fenti problémát természetesen csak valós változókat tartalmazó alakban is meg-—

fogalmazhatjuk, ilyenkor a [23/ feladat második feltételrendszerét a vele egyenértékű

din—d,) : 0 (i : n—l—1, . . ., N) /24/'

feltételekkel is helyettesíthetjük.

Az optimális r számú mérési helyet pedig azon y,- pontok adják, amelyekre d, : 1,, valamint az optimális célfüggvényérték szolgáltatja a hozzá tartozó variancia—értéket.

3. Adott kutatási biztonság optimális elérése

Mint azt már a bevezetőben is említettük, természetesen vetődik fel az előbbiekben vázolt feladat megfordítása is: előírt korlátnál kisebb variancia-érték biztosítása minél kevesebb újabb mérési hely figyelembevételével.

Ennek a feladatnak a megoldását, azaz a minimális r értéket és a hozzá tartozó mérési helyeket megadó módszert is a /23/—hoz hasonló vegyes programozási feladattal adhatjuk

meg. Ez a feladat /23/—tól csak abban különbözik, hogy a célfüggvény értékére felső korlátot

adunk, ez felel meg a variancia-értékre adott korlátozó feltételnek, valamint az új célfügg- vény az r-et adó harmadik feltétel bal oldala lesz, hiszen r adja a további mérési pontok számát, amelyet minimalizálunk. Tehát feladatunkat a

dfzi (i:1,....n) d,E(0,1] (i:n-l—1,...,N)

N

j;dj/ljz1

N

diu'l'j2yijdídjajz Jim (í:1'Z!---1N)

N

m— j); 1; el,-no § (fd-7700 /25/*

vegyes egész valós programozási probléma alakjában is megfogalmazhatjuk, ahol ismét y,, 1,- és d,- (i : 1, . . ., N) az ismeretlenek. A programozási feladat diszkrét jellege most is megszüntethető azáltal, hogy a második feltételrendszert a /24/ feltételekkel helyettesít—- jük, és a di változókat is valós, folytonos változóknak tekintjük.

4. Az eljárás numerikus bemutatása

A /23/ programozási feladat konkrét alkalmazására a következő oldalon levő táblában!

alappontokkal és értékekkel megadott függvényt választottuk.

Az egyszerűbb számolások érdekében az r : 1 esetet, azaz egyetlen újabb mérési hely' kijelölését vizsgáltuk. A megadott alappontokon kívül egy ekvidisztáns ponthalmazt válasz——

tottunk, mint a lehetséges újabb mérési helyeket. Ezek a pontok a következők voltak:

1 , 1 1 . 1 .

(í 'l'l'ír í'l']"í)v l26/

(6)

KlSÉRLETEK TERVEZESE 1 1 19

ahol i, ] egymástól függetlenül befutották a (0, 1, 2, 3, 4, 5] indexhalmazt. A y(h) variogram-

ként a szferikus formát választottuk.

A függvény alappontjaí és értékei

l

X2 0 1 /3 2/3 1

X1

o 0 0 0 o

1 /3 0 0.110 883 0220 398 0327195 2/3 0 0220 398 0429 956 0,618 370 1 0 o,327195 0,618 370 0841 471

Esetünkben tehát n : 16, N : 16—936 : 52, így a /23/ feladat összesen egy ,a, valamint 52 Á,— és 36 d,- ismeretlent tartalmaz, azaz az ismeretlenek együttes száma nagy:

14—524—36 : 89. Minthogy i: 1, . . ., 16 esetén d,- : 1, az ismeretlenek számát 16-tal

csökkentettük oly módon, hogy az első tizenhat da értékét, az 1—eket a modellbe beírjuk.

Ekkor a modell a következőképpen módosul:

dja—dj) : 0 (j : 17, 52)

52 16

2 MH" Z A; :1

;:17 j—1

52

d,. : 1 /z7/

j:17

16 52

H—l— Z yüAj-l— Z yi, djlj : 710 (i : 1, . . ., 16)

1-1 1-17

18 52

díy-l— 121 yi; dalia—127 y,.j d,d,A,——d,.7w : 0 (i : 17, ..., 52)

a 'n]

16 52

PH" Z Áj7jo'l' Z A; (Wie—* mi"

fal j-17

hiszen a célfüggvényben a 7700 tag konstans, így elhagyható. Meg kell viszont jegyeznünk, hogy az így felirt célfüggvény optimális értékéből yoo-át még le kell vonnunk, hogy az opti—

mális variancia-értéket megkaphassuk.

A [27/ programozási probléma az így felírt alakban nem lineáris, folytonos változókkal rendelkező feladat. A célfüggvény kvadratikus, viszont az ötödik feltételrendszer harmadik tagjában szereplő összeadandók harmadfokúak, hiszen itt di, dj, Áj egyaránt ismeretlenek.

így a feladatot általános alakú, nem lineáris, nem konvex problémákra kidolgozott mód——

szerekkel oldhatiuk meg. A konkrét megoldást az Arizonai Egyetem Számítóközpontiának Cyber 175 típusú gépén végeztük el. A teljes számítás 352 másod perc gépidőt vett igénybe.

Az újabb mérési hely figyelembevétele nélkül a függvény átlagértékeként O,195 788- adódott, ezen becsült érték varianciájára pedig a (19) összefüggés alapján 0,001 258. Az opti—

mális megoldásban egyedül a (0,25, 0,25) ponthoz tartozó d,- érték adódott 1 értékűnek, az összes többi dj érték pedig O volt, azaz az optimálisan megválasztott újabb mérési hely a (0,25, 0,25) pont. Az optimális célfüggvényértékből 7700 értékének levonása után megkap—

tuk az újabb pont figyelembevételével adódó csökkentett variancia-értéket is: 0,001 162,

azaz egyetlen újabb pont és mérés felvételével a becslés bizonytalansága legkedvezőbb

esetben 0,000 O96—dal csökkenthető, és a legkedvezőbb csökkenés a (0,25, 0,25) pont felvételével és az ezen pontban történő méréssel érhető el.

(7)

1 120 DR. SZIDAROVSZY FERENC

5. A modell továbbfejlesztési lehetőségei

A továbbiakban néhány olyan kérdést említünk meg, amelyek az itt bemutatott modell továbbfejlesztésével kapcsolatosak.

A) Az univerzális Kriging és a Cokriging módszer esetére minden további nélkül átirható az eljárás. A részletek mellőzésével megjegyezzük viszont, hogy a Cokriging mód—

szer esetében (10) egyszerre több paraméter becslését végezzük el. így a 2. pontban leirt optimumfeladat többcélúvá válik, mivel az egyes célokat a különböző paraméterek átlag—

értékeinek megbízhatóságát jellemző varianciák adják. Ekkor valamilyen többcélú progra- mozási algoritmust kell a megoldáshoz választanunk. A 3. pontbeli problémánál pedig mindegyik varianciára korlátot kell adnunk.

B) Főleg geostatisztikai kutatások esetén gyakran előfordul, hogy nem egy-egy para- méter átlagértékét kívánjuk becsülni, hanem a paraméterek értékeit egy előirt ponthal—

zmazon. Például ez a feladat merül fel a szintvonalak meghatározásakor. llyenkor a célunk nem a paraméterek átlagértékei szórásának minimalizálása, hanem az előirt alappont—

halmazon felvett értékek szórásai közül a legnagyobbnak a minimalizálása. Ha Dl, . . ., D, jelöli az egyes pontbeli szórásnégyzeteket, akkor a legnagyobbnak az értékét E-vel jelölve

nyilvánvalóan a

D6§E (i:1,...,r) további feltételek mellett az

5—3 min

szélsőérték-feladatot kell megoldanunk. Ha több paraméter együttes figyelésére is alkal—

massá akarjuk tenni a modellt, akkor minden egyes paraméterhez egy-egy EU) értéket vezetünk be az előző módon, és ismét többcélú programozási módszerrel kísérletezünk.

C) A bemutatott modell minden további nélkül átírható oly módon, hogy még a következőkre is tekintettel legyen. Ha ugyanis az előző esetben egy-egy pontbeli para- méterérték becsléséhez nem kívánjuk összes mérési; értékeinket feihasználni (például a variogram szerint a tőlük távoli értékek már függetleneknek tekinthetők), akkor a pont- helytől függően is megállapodhatunk bizonyos d,- értékek zérus voltában. Ezzel tudjuk ugyanis garantálni, hogy a neki megfelelő pont a feldolgozásból kimaradjon.

D) A (25) és a (23) feladat vegyes programozási problémát jelent, viszont a di(1——-d,-) : 0 feltételekkel folytonos feladatként kezelhetők. Természetesen ugyanez elmondható az univerzális Kriging és a Cokriging módszer esetében is. Viszont észrevéve, hogy a varian- ciáknak megfelelő képletekben a ti,-dili harmadfokú tagok szerepelnek, egyszerüen, bár afeladat dimenziójának megnövelésével, kvadratikus programozási feladatra redukálhatjuk.

Ha ugyanis 3

Díj : did] (Vw)

változót is ismeretlenként kezeljük, ez az újabb feltétel is másodfokú, viszont a dídej : Díj-Aj

egyenlőség következtében a programozási feladat harmadfokú tagjai is másodfokúakká

válnak.

jelen dolgozatban költséges mérések optimális tervezésére mutattunk be nem lineáris programozási feladatra vezető módszert. Az általános modell adott számú mérési pont

(8)

KlSÉRLETEK TERVEZESE

1 121

optimális tervezésére, valamint adott biztonsági szint eléréséhez minimálisan szükséges mérési pontok számának és maguknak a pontoknak a meghatározására alkalmas. A javasolt eljárás alkalmazhatóságát konkrét számpéldán is bemutattuk. Az eljárás továbbfejlesztési lehetőségei között megvizsgáltuk például a Cokriging és az univerzális Kriging módszert, valamint izovonalak szerkesztéséhez szükséges optimális fúrási stratégia kialakitásának módját. A módszer matematikailag vegyes programozási feladat megoldását igényli, amelyet kvadratikus programozási problémára vagy harmadfokú, kisebb dimenziójú feladatra is redukálhatunk.

IRODALOM

251 2 (Dlgelhomme, ]. P.: Kriging in the hydrosciences. Advances in Water Resources. 1978. évi 5. sz.

-— 66. o .

(2) Gomba/ati, G. — Volpi, G.: Groundwater contour mapping in Venice by stochastic interpolators.

1. Theory. Water Resources Research. 1979. évi 2. sz. 281—290. old.

(3) Matheron, G.: The intrinsic random functions and their applications. Advance: in Applied Pro- babilíty. 1973. évi 5. sz. 439—468. old.

(4) Szép Jenő: Analízis. Matematikai ismeretek gazdasági szakemberek számára. Közgazdasági és Jogi Könyvkiadó. Budapest. 1972. 777 old.

(5) Krekó Béla: Lineáris algebra. Közgazdasági és Jogi Könyvkiadó. Budapest. 1976. 573 old.

(6) Szidarovszky Ferenc: Bevezetés a numerikus mádszerekbe. Közgazdasági és Jogi Könyvkiadó.

Budapest. 1974, 387 old.

(7) Szidarovszky F. — Yakowitz, S.: Principles and procedures of numerícal analysis. Plenum Press.

New York. 1978. 331 old.

(8) líorgó Ferenc: Nemkonvex és diszkrét programozás. Közgazdasági és Jogi Könyvkiadó. Budapest.

1978. 436 0 d.

(9) Optimumszámitási modellek. Főszerk.: Kósa András. Műszaki Könyvkiadó. Budapest. 1779. 865 old.

(10) ]ournel, A. G. -— Huíjbregts, Ch. J.: Mining geostatistícs. Academic Press. London —- New York.

1978. 600 old.

PE3IOME

ABTop geMoncrpnpye'r meron, npnroAHbli—i Ana pemenua aagau Henunem—ioro npurpum—

anoaam—m, cansanubix c OHTHMaanbIM l'll'löHMpOBöHHeM Aoporocronumx oőcnenoaam—iü. 06- 1145; mortem, npuroAHa Ann onpeneneuun ManMaanoro Konuuecraa M3MepHTeanle Touek

" caMnx roueK, Heoőxonwmmx nna nonyuemm AaHHoro uucna namepurenbnux rouek. Aarop M Ha KOHerTHOM npumepe nonaabiaaer anMeHHMOCTb npegnaraemoro Merci.-ia. CpeAn Bos- Momnocreü nanbueümero paasmm aerop npuaonur Me'ron prrmira " Kakpnmnra, a rak—

me cnocoő oőpaaoaanua onmmanbnoü Apeneaoü crparernn, Heoöxonumoü nm: nocrpoem—m usanuHuí—i.

Mei-op, paccuwran Ha pewei-me Maremaruuecm cmetuam-loü sanaun nporpaMMupoaaHm, Koropyio Momno csecm Takme H Ha ksagpamuecxyio, Kyőw-iecxyio sanauy nporpaMMupo- senna menbmnx paamepoa.

SUMMARY

The study presents a method for solving non linear programming tasks aimed at optimal planning of expensive surveys. The general model can be used for determining the minimum number of measuring points needed for obtaining a given number of measuring points as well as for determining the points themselves. The author shows the applicability of the procedure put forward with factual example. From among the possible improvements of the procedure the Kríging and Cokriging methods and the elaboration of optimal boring strategy for drawing isolines a_re given.

From mathematical point of view, the method reauires the slaving of a mixed pro—

gramming task which can be reduced to auadratic or cubic ones.

.5 Statisztikai Szemle