NOT FOR QUOTATION WITHOUT THE PERMISSION OF THE AUTHOR
THE EXPEXIIENTAL DESSGN
OF
AN OBSERVATION IKETWORK:SOFTWARE
AND
EXAMPLESK Fedorov, S. Leonov, M. Antonovski, S. P i t o v r a n o v
January 1987 WP-07-05
Working R a p e r s are interim reports on work of t h e International Institute f o r Applied Systems Analysis and have received only limited review. Views or opinions e x p r e s s e d h e r e i n d o not necessarily r e p r e s e n t t h o s e of t h e Institute or of i t s National Member Organizations.
1NTE;RNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS 2361 Laxenburg, Austria
Foreword
One i m p o r t a n t theme within t h e I I A S A Environmental P r o ~ r a m i s t h e optimiza- tion of monitoring systems. This Working Paper i s a contribution to t h a t activity.
Dr. Leonov was a n a s s o c i a t e r e s e a r c h s c h o l a r at I I A S A from November 1986 t h r o u g h J a n u a r y 1987. His a r r i v a l t r i g g e r e d some p a r t i c u l a r l y f r u i t f u l discussions o n t h e design of monitoring systems, and led to t h e development of t h e system re- p o r t e d in t h i s p a p e r . Although t h e example used to test t h e a p p r o a c h is r a t h e r e l e m e n t a r y f r o m t h e standpoint of a pollution c o n t r o l a g e n c y , t h e i d e a s p r e s e n t e d p r o v i d e a s p r i n g b o a r d f o r f u r t h e r s t u d i e s at I I A S A a n d in Dr. Leonov's home insti- tution in Moscow.
Acknowledgement
W e are most ~ r a t e f u l to Professor R.E. Munn f o r his constructive criticism and recommendations which were extremely helpful in the preparation of this pa- per.
THE
EWQUKEHTAL DESIGNOF
AN OBSERVATIONAL NETWORK:SOFTWARE
AND
EXAMPLESV Fedorov, S. Leonov, M . Antonovski, S. Pitovranov
I. INTRODUCTION
Because of i n c r e a s i n ~ anthropogenic impacts on t h e biosphere, t h e r e h a s been a r a p i d development of environmental pollution monitoring systems. Durlng t h e l a s t 15-20 y e a r s many monitoring programs h a v e been initiated on t h e national (including c o u n t r i e s such as Canada, G r e a t Britain, USA, USSR and o t h e r s ) and international levels (including t h e European Monitoring and Evaluation Program (EMEP) under t h e auspices of t h e European Commission f o r Europe, and t h e Back- ground Air Pollution monitor in^ Network under t h e auspices of t h e WMO). Environ- mental monitoring systems have been defined in s e v e r a l ways. An Intergovernmen- tal Working Group defined monitoring as "a system of continued observation, meas- urement and evaluation f o r defined purposes" (International Working Group 1971).
Yu. Izrael (1974) h a s provided t h e following definition:
-
observations of t h e state of t h e environment and f a c t o r s affecting t h e en- vironment;-
assessment of t h e state of t h e environment and impact f a c t o r s ;-
prediction and assessment of t h e f u t u r e state of t h e environment.According t o Izrael, observations are only t h e f i r s t s t a g e of a monitoring system.
Anthropogenic impacts on t h e b i o s p h e r e h a v e different spatial and temporal scales. That i s why monitoring systems also have to b e classified by d i f f e r e n t spa- tial levels, h e r e w e define t h e Levels as local, regional and global.
Local monitoring systems a r e usually r e q u i r e d in u r b a n a r e a s (city, industrial complex) to d e t e c t anthropogenic impacts o n human health and t h e environment.
The geographic s c a l e i s of t h e o r d e r of s e v e r a l kilometers to one hundred kilome- t e r s . The monitoring systems assess environmental quality, which is usually based on such c r i t e r i a as maximum permissible pollutant concentrations in a i r , drinklng water, e t c . Also t h e y provide d a t a f o r t h e assessment of economic losses resulting from environmental pollution. Usually, i t i s possible t o define relationships between emissions and concentrations of pollutants o n local levels.
Regional monitoring systems s e r v e a similar purpose, but o p e r a t e o v e r a con- siderably l a r g e r s c a l e (approximately 1000 km). A major problem i s t h a t s o u r c e - r e c e p t o r pollution relationships are much more complex on t h e regional level than they are on t h e local level.
The goal of global monitoring systems i s t o o b s e r v e and a s s e s s changes in t h e biosphere on a planetary or hemispheric s c a l e . In t h i s c a s e anthropogenic im- p a c t s a r e o b s e r v e d as averages o v e r t h e e n t i r e world community. Accordingly, s o u r c e - r e c e p t o r linkages a r e even more complicated than at t h e regional level.
The impacts of anthropogenic s o u r c e s are mediated by p r o c e s s e s having t h e fol- lowing specific f e a t u r e s (Rovinski and Wirsma, 1986). First, such p r o c e s s e s a r e
a s s o d a t e d with pollutants having small ooncentrations which =use large-scale cu- mulative effects. Second, these processes are always associated with long-range transport of pollutants, mainly in t h e atmosphere, but sometimes in t h e hydro- s p h e r e . Third, global impacts from anthropogenic activities have long lead times before they are noticed. Therefore, i t seems t h a t even in f u t u r e studies, global monitoring programs will focus on statistical analysis of observation data (see, f o r example, Izrael, Rovinski, Antonovski et al. 1985), and not on defining souroe-receptor relationships.
According to t h e above-mentioned definition of Izrael, t h e objectives of moni- toring include prognoses of environmental states. The tool for such forecasts i s mathematical models of different kinds (from p a r t i a l differential equations to sim- p l e regression equations). There are many obstacles to t h e application of mathematical models in environmental sciences. W e mention h e r e only a few of them: t h e lack of knowledge of processes and physical parameters of ecosystems;
t h e complexity of physical processes in environmental systems; t h e stochastic c h a r a c t e r of various processes (for example, atmospheric turbulence); measure- ments of different environmental parameters with Large e r r o r s ; etc.
The goal of this study i s to suggest a p r o c e d u r e which can be useful for optim- izing t h e design of monitoring networks. The approach and t h e numerical algo- rithms presented in t h i s p a p e r are to s o m e e x t e n t a continuation of t h a t by Fedorov, 1986. Here, however, emphasis will be placed on t h e technical aspects of t h e numerical algorithms relating to the construction of optimal observational net- works.
The proposed versions of algorithms are oriented to solving problems of t h e optimal location of observation stations In a given region. Their generalization to o t h e r optimal experimental design problems should b e straightforward. Any princi- pal changes to b e made would b e in the block describing "controllable area." As anticipated r e a d e r s will be from t h e environmental pollution or meteorological paradigms, all illuminative examples are r e l a t e d to a i r pollution monitoring.
Wewish to emphasize two main assumptions which are crucial to the approach.
1. The optimal design of an observational network is model oriented. I t is as- sumed t h a t t h e t r e n d s of t h e observed values oan be (at least approximately) described by a mathematical model containing unknown parameters. For in- stance, it could be a Gaussian plume model where diffusion coefficients, a n af- fective s t a c k height and emission intensity have to be estimated (S. Hanna et d., 1982).
2. All uncertainties (observational errors, fluctuations of processes under in- vestigation, small imegularities, deviations of t h e model from t h e "true"
behavior, etc.) are absorbed by additive e r r o r s , which are assumed to b e random.
Assumptions about the randomness of errors are crucial to t h e whole ap- proach because all objective functions (both in analysis and design) are for- mulated as expectations of some deviations of estimators from t h e "true"
values. Most frequently i t is the variance of a n estimator or the variance- covariance matrix and s o m e functions of i t in multidimension cases.
A summary of 1 and 2 leads us to the following presentation:
Yi
=
q ( z f ~3) + ti 9 (1)where yi is an observed value at a n 1-th station; zf ere t h e coordinates of t h i s sta- tion; q ( z ,9) is an a p r t o f i given function (for instance, i t could b e a concentration of some pollutant); cf is a random d u e with z e r o mean ( E ( r i )
=
0). Usually, q ( z ,iP) is calledthe
"response function."The algorithms presented in this paper are oriented to the case where e r r o r s of observations are uncorrelated: E [ci c,]
=
bZA"(zi)bi,, where A(z) is t h e so- called "effectiveness function" reflecting t h e accuracy of observations at t h e given point z.
I t is assumed throughout this paper t h a t yi is a scalar. The generalization f o r more complicated situations, f o r instance yi e i t h e r a vector o r a function of time, is straightforward (compare with Fedorov, 1972, Ch.5; Mehra and Lainiotis, 1976).
One can apply t h e method to a vector case when t h e concentration of several pollutants have t o be observed. If the dynamics of some environmental charac- teristics are of interest then i t becomes necessary to consider responses belong- ing t o some functional space.
IL
0- CRITERIAThe algorithms presented in this paper comprise two main types of optimality criteria: t h e f i r s t is related to the variance-covariance matrix of estimated parameters, while t h e second is based on variance characteristics of t h e response function estimators. Details can b e found in Fedorov, 1972; Silvey, 1980; Atkinson
& Fedorov (to b e published).
Table 1 contains optimality c r i t e r i a which can b e handled with the help of t h e software described later. The c r i t e r i a used in t h e statistical literature are in t h e f i r s t column of t h e table; formal definitions of optimality c r i t e r i a are in t h e second; and t h e corresponding dual optimization c r i t e r i a are formulated in t h e third column (see, f o r instance, Atkinson & Fedorov (to b e published)).
TABLE 1.
optimarity c r i t e r i a
D-criterion generalized D-criterion
linear criterion
A(zIfT(= )023 (2)-tr D , A ( z ) ~ ~ ( z ) r n f ( z )
-
f r AD,can b e transformed to t h e previous case with
extrapolation
I
ci ( l o .0
I I
A ( z ) ~ ~ ( z ) D ~ ( z ~ ) I ~ ~ ( z ~ , ~ ) ~1
~ = f ( z ~ ) f ~ ( z ~ ) .Theoretically all of t h e algorithms discussed are valid for t h e case of Linear parametrization:
where j ( z ) i s a v e c t o r of given functions. How to handle nonlinear models will b e considered in Example 3.
The following notations are used in Table 1:
-
D=
D ( O=
N D ( ~ ) , where D=
D ( O i s a normalized v a r i a n c e - c o v a r i a n y matrix, D(*) i s a variance-covariance matrix of t h e least s q u a r e estimator 9 ,n
D - l ( t )
=
M(4)=
) f ( z ) f '(21
or=
pi A1"(zi (zi1.
X i =l
- 4
i s a design, i.e.,t=
fpi ,zi , where pi i s a fraction of observations which h a s to b e located at a point zi; pi could be t h e duration, frequency or t h e precision of observation;-
m i s a number of unknown p a r a m e t e r s (dimension of 9 ) ;-
d (2.4)=
j (Z )Dj ( z ) is a normalized variance of t h e estimator q ( z ,s)
at agiven point x ;
-
X is a controllable region, zi E X ;-
A i s a utility matrix, usually reflecting t h e significance of some p a r a m e t e r s or t h e i r l i n e a r combinations;-
o ( z ) is a utility function, usually reflecting t h e i n t e r e s t of a p r a c t i t i o n e r in t h e value of t h e response function at a point z.The existence of a nonsingular optimal design is assumed f o r all optimality c r i t e r i a in Table 1. Singular optimal designs (i.e. a n information matrix ~ ( 4 ' ) i s singular, IM(4') ( = 0, in t h e r e g u l a r case D(4) = M-l(t)) can o c c u r when rank A
<
m.
In p r a c t i c e one c a n easily avoid singular designs applying to t h e regular- ized version of t h e initial problem (see Fedorov, 1986, section 2):Q,[D(t)l
=
Q [ I ( ~ - P ) M ( ~ )+
pM(tO)ll l .
(3)where
I
M(to)I +
0.Objective function (3) c a n also b e used in cases where it is necessary to com- plement existing networks defined by
4,
by some new observational stations. D- and A-criteria are usually used when all unknown parameters are equally of i n t e r e s t . The f i r s t o n e i s preferable, being invariant to linear transformation of unknown parameters (for instance when o n e needs to rescale some of them). This assertion c a n b e easily verified with t h e help of t h e corresponding dual optimization prob- lems (Table 1, column 3). When some parameters are more important than o t h e r s , matrix A i s usually chosen diagonal with elements Aa,(a ;=
1 ) reflecting t h esignificance of the corresponding parameters $,(a
= -
1,m).The last two c r i t e r i a can b e used when a n experimenter is interested in t h e explicit estimation of a response function q ( z , 6 ) . For instance, if t h e r e are points z1,z2, ...,z, of special interest, then o ( z )
= E
d(z 7 k ) , where d ( z 7 ' ) i s 6-3 t = l
function, and O(D)
=
C d (zt ,€I.t = l
IIL FIRST-OWER I'rEBATlVE A L G C ) ~ 111-1. ?he a l g o r i t h m
W e start with t h e iterative algorithm of t h e following form ( f o r details see Fedorov 1986):
where
4,
is a c u r r e n t design on a s t e p s ,ts =
Izis ,pi,, i = l , n s- 1,
C pis=
1 ,i =I
& =
tzar; i= -
l , n s1
is a supporting set of t h e design;t(zs
) i s a design with t h e measure totally located at a point z,.
For designs with a finite number of supporting points, formula (4) means that P(,,+l
=
( l - , ) ~ ~ ~ , i d * , Pi.,s+l=
(1-a,I~~.,+a,, HZ, =zt,.,'The algorithm provides so-called forward and backward procedures. In the back- ward procedure, the 'least informative" points are deleted from t h e c u r r e n t design, while conversely t h e forward procedure includes t h e new, "most informa- tive " ones.
111-2. Selection of Izs
1
andI
a,1.
For t h e forward procedure:
For t h e backward procedure:
p,'
=
p (zS3 i s a weight f o r a point z,The algorithm provides t h r e e choices of gain sequence lysj:
-
1(a) 7,
- -
, s =1,2,...
; no i s a number of supporting points in a n initial n 0 + sdesign. With t h i s choice of 7,. one c a n simulate t h e subsequent inclusion (deletion) of t h e m o s t (least) informative stations.
(b) 7, i s defined by t h e s t e e p e s t descent method, which provides t h e l a r g e s t d e c r e a s e of t h e objective functions in t h e chosen direction t ( z ).
(c) 7,
=
C o , where Co is a small constant (0.01+
0.1) which i s defined by a u s e r . This sequence does not satisfy traditional conditionswhich are usually implied to p r o v e t h e convergence of t h e i t e r a t i v e algorithms, but may b e useful f o r t h e construction of t h e d i s c r e t e designs.
Numbers of steps (length of excursion) f o r t h e forward and backward p r e c e d u r e s (Mar and n b a c respectively) are defined by t h e user.
111-3. D c r i t e r i o n .
The algorithm '?)OPTJ' i s oriented f o r t h e construction of D-optimal designs providing t h e minimum of t h e determinant ( D ( 9 ) ( , D ( 9 ) i s a covariance matrix of t h e parameters' estimators. Geometrically t h e minimization of t h e determinant i s equivalent to t h e minimization of t h e volume of t h e ellipsoid of concentration f o r t h e parameters' estimators. Simultaneously t h e algorithm minimizes sup X(z )d ( z , t ) ( s e e Table 1 ) securing a n effective estimation of t h e response f unc-
r EX
tion o v e r set X. Moreover, in t h e case of normally distributed e r r o r s ci D-optimal design e n s u r e s t h e b e s t value of t h e noncentrality parameter when t h e hypothesis
s u & ~ ~ ( z . * , )
+
d , d>
0 .i s tested (see Fedorov, 1986).
A function q ( z , t s ) h a s t h e following presentation f o r t h e D-criterion (see Table 1):
where d ( x , # , )
=
f r ( x ) ~ - l ( # s ) f ( x ) l s e e a l s o page 4.The formulae f o r sequential recomputation of t h e c o v a r i a n c e matrix and t h e determinant a r e
111-4. Some n o t e s o n t h e d g o t i t h m . S t o p p i n g r u l e
The calculations are terminated if
(a) t h e convergence c r i t e r i o n i s attained f o r t h e forward procedure:
I
cp(z,+)I
<
6, w h e r e 6 i s defined by a u s e r (this means t h a t t h e value of t h e mdirectional d e r i v a t i v e is s m a l l e n o u ~ h and, subsequently,
ts
i s close e n o u ~ h t o t h e optimal design).(b) a given number of i t e r a t i o n s i s attained.
Merging of s u p p o r t i n g p o i n t s i n t h e firward p r o c e d u r e
Let hk b e a size of t h e k-th grid element defined during t h e mapping of X, k
=
; L is a dimension of controllable region X. Ift h e n a point xi is merged with a point x i , constant C,,, being defined by a user.
Deleting of p o i n t s with s m d l w e i g h t s i n t h e f o r w a r d procedure.
If f o r some i , pis,
<
6 t h e n a point xi,, is deleted from t h e design. The covariance matrix a n d c r i t e r i o n value are recomputed, formulae (5), (6) being used with a= vi ,,
a n dBoth l a t t e r p r o c e d u r e s h e l p t o avoid designs with a Large number of supporting points.
lV.
OPTIMIZATION ALGORITHM OFTHE EXCHANGE TYPE
The algorithm h a s t h e form
ts
=ts
+ as t( z s ) where a, can b e e i t h e r positive o r n e ~ a t i v e .From a computational point of view, t h e main d i f f e r e n c e in algorithm ( 7 ) from t h e one described in Section 3 i s t h a t t h e whole design i s not recomputed at e a c h s t e p ; all modifications c o n c e r n only newly included (a,
>
0) or deleted (a,<
0) points, which explains t h e origin of t h e t e r m "exchange" in t h e t i t l e of t h e algo- rithm (see a l s o Fedorov, 1986). The various modifications of t h e "exchange type"algorithm are p a r t i c u l a r l y useful when some subset of a n initial design h a s to b e included in the final design (some p r e s c r i b e d observational s t a t i o n s h a v e to b e included in t h e final o b s e r v a t i o n a l network). The algorithm c a n b e easily a d a p t e d to solve t h e regularized v e r s i o n s of t h e originally singular design problems con- s e r v i n g some "regular" f r a c t i o n s of a n initial design.
The presented s o f t w a r e contains three modifications of t h e exchange pro- c e d u r e .
(1) Deleting t h e l e a s t i n f o r m a t i v e p o i n t s f r o m t h e initial d e s i g n .
The backward p r o c e d u r e is executed (some points are deleted) with a,
=
-I/ n o , n o i s t h e number of points in t h e initial designA number of s t e p s f o r deleting ( m i t e r
=
nbac) i s chosen by a u s e r . All points in t h e final design (normalized to: Cpi=
1 ) h a v e equal weights:pi
=
I / (no--nbac), n b a c<
n o .This p r o c e d u r e c a n b e used, f o r instance, when i t i s necessary to find a n d remove a given number of t h e l e a s t informative s t a t i o n s (see example l ( b ) .
(2) I n c l u s i o n of t h e most i n f o r m a t i v e p o i n t s . The forward p r o c e d u r e i s executed with
z,
=
zf=
Argmin q ( z ,#,
);ta
a number of s t e p s f o r including (miter
=
nfor) i s chosen by a u s e r ; final weights are equal toF o r both of t h e a b o v e p r o c e d u r e s , t h e normalization of t h e c o v a r i a n c e matrix i s c a r r i e d o u t during t h e last s t e p , as otherwise C p i ,,
<
1 f o r t h e backward pro-i
c e d u r e , and C pi ,,
>
1 f o r t h e forward procedure.i
Normalization i s not e x e c u t e d during t h e intermediate s t e p s in o r d e r to make tan- gible e i t h e r t h e d e c r e a s e of t h e determinant ID(#,)( due to t h e deletion of t h e observational stgtions or i t s i n c r e a s e due to t h e inclusion of stations.
(3) Standard e z c h a n g e p r o c e d u r e .
Forward and backward p r o c e d u r e s are executed subsequentially, t h e initial p r o c e d u r e being chosen by a u s e r . H e r e , t h e number of s t e p s f o r t h e forward and backward p r o c e d u r e s are equal:
W r
=
nbac=
nn,and t h e maximal number of i t e r a t i o n s h a s t h e form miter
=
2 n n *ko, k o i s a n integer.A measure [ i s a probability measure at t h e end of t h e ' l a r g e iteration", t h e length of which is equal to 2 n n ; t h i s f a c t explains t h e choice of a value m i t e r .
The c h o i c e of
12,
j i s as d e s c r i b e d above,ys forward p r o c e d u r e -in (7, , p ; ) , back ward p r o c e d u r e T h e r e are t w o v a r i a n t s f o r t h e choice of gain sequence 17, j:
-
1's
-
no+l+l s = 1 2 ; 1 i s a n i n t e g e r p a r t of ( s -I)/ 2 n n ;ys c h a n g e s a f t e r executing both forward and backward p r o c e d u r e s , i.e., i t i s a ' l a r g e iteration";
(b) ys
=
Co ,Co i s a c o n s t a n t defined by t h e user.The popular Mitchell algorithm (Mitchell, 1974) can b e considered as a p a r t i c u l a r case of t h i s version. I t is w e l l known t h a t t h e Mitchell algorithm does not generally converge to a n optimal solution.
Some n o t e s o n t h e a l g o r i t h m . Stopping r u l e
( a ) The c o n v e r g e n c e c r i t e r i o n is attained at t h e last s t e p of t h e forward prc+
c e d u r e
(b) t h e maximal number of i t e r a t i o n s i s attained.
Merging of p o i n t s for f i r w a r d procedure This i s t h e s a m e as d e s c r i b e d in Section 111.
n z i n g of initial p o i n t s .
Some points in a n initial d e s i ~ n may b e fixed by a u s e r (locations and weights).
E n s u r i n g that the w e i g h t s a r e positive.
If for t h e backward p r o c e d u r e
<
ys f o r some s (see (8)), then f o r t h e sub- sequent forward p r o c e d u r ea,,
=
p i , f o r some s'.
This modification i s n e c e s s a r y to k e e p
tS
in t h e set of probability measures at t h e end of t h e ' l a r g e iteration".V.
UNEAROPTIMALITY CBITEBUAlgorithms LINOPT and LINEX a r e intended f o r t h e construction of l i n e a r optimal designs providing minima of t h e value
where A i s a utility matrix (it i s a symmetric nonnegative definite ( m Xm ) matrix) chosen by t h e u s e r according to his needs.
The major difference in t h e algorithms LINOPT ( f i r s t ~ r d e r i t e r a t i v e algo- rithm) a n d LINEX (optimization algorithm of t h e exchange type) from DOPT a n d DOPTEX respectively, i s t h a t t h e function ~ p ( z ,[, ) h a s t h e following presentation:
VI. USER'S W I D E
1 . Mapping o j a controllable r e g i o n X Program MAP i s intended f o r mapping of a controllable region X. The c u r r e n t version of t h e program handles one- and two- dimensional regions but generalization to l a r g e r dimensions is not difficult. The region X i s defined on a uniform g r i d with given densities f o r e a c h variable. Such a presentation of X is explained by t h e f a c t t h a t usually a u s e r deals with i r r e g u l a r regions which cannot be d e s c r i b e d analytically (non-convex, with subregions where t h e location of observational stations is impossible, f o r example, lakes, densely populated a r e a s , etc.). Two output files are c r e a t e d by t h e program:
"REG.DATN contains t h e d a t a in i t s original scale,
'SCALE.DATn contains t h e normalized d a t a (-1
s zInor)
S +1, i = l J ).2. Programs DOPT, DOPTEX, LINOPT, LINEX, utilize t h r e e files:
"0UT.D AT" i s f o r output information (see examples),
"REG-DAT" contains t h e design's grid,
"DES.DATn contains a n initial design.
3. The s t r u c t u r e of a v e c t o r of basic functions f ( z ) must b e set i n t h e subroutine resp:
where m i s t h e number of unknown p a r a m e t e r s ,
I =
V1(2),. . -
# f , ( ~ ) ) ~ . . z= .
P.
If t h e effectiveness function X(z) is not constant, then instead of j , ( z ) t h e func-
-
tion ~ " ~ ( z ) f , ( z ) h a s t o be programmed, a
=
1,m.
4 . All auxiliary subroutines (matrix inversion, calculation of t h e initial d e t e r - minant, minimization of a function ~ ( z ,[, ) e t c . ) f o r programs DOPT.DOPTEX and LINOPT.LINEX a r e saved in t h e files 'SUBD .F" and 'SUBL.Fn, respectively.
VII. Examples
E z a m p l e 1. L i n e a r p a r a m e t t i z a t i o n , D c t i t e t i o n .
To illustrate t h e possibilities of t h e proposed software, l e t us consider a com- paratively simple example based on a i r pollution d a t a from Modak and Lohani, 1985. The p a r t i c u l a r example we shall use is shown in Figure 1, which gives iso- p l e t h s of monthly mean values of SOz concentration f o r 9am in Taipei City, Taiwan.
The original network contains eleven observing stations (see Figure 2). To begin.
t h e underlying model w a s chosen as a polynomial of t h e second d e g r e e with u n c o r r e l a t e d random additive e r r o r s :
where ( z l f , z z i ) are c o o r d i n a t e s of t h e i-th station. Of c o u r s e , t h i s model i s t o o simple f o r a good approximation of t h e p a t t e r n presented in Figure 1 , but because of i t s simplicity one can easily understand t h e main f e a t u r e s of t h e software.
The optimality c r i t e r i o n w a s taken equal t o t h e normalized determinant of variance-covariance matrix (D-criterion).
(a) Completely new network. The purpose of t h i s algorithm i s t o find t h e
"best observation" network under t h e assumption t h a t t h e r e are no constraints on t h e number of stations and t h e i r locations e x c e p t t h a t t h e stations have to b e within t h e city's a r e a .
The r a t i o of determinants f o r t h e original and optimal locations i s g r e a t e r than 1 0 (see Printout 1 ) . One can o b s e r v e (Figure 3) a typical ( f o r t h e conven- 4 tional optimal design) location of observation stations: m o s t of them have to be on t h e boundary of t h e a r e a and only a few (in o u r c a s e only one) inside it. This should be compared with t h e r e s u l t by Modak and Lohani, 1985, p.14, based on t h e so-called "minimum spanning t r e e " algorithm, where observing s t a t i o n s a r e mainly located inside t h e area.
However, a comparison of r e s u l t s is conditional since t h e a u t h o r s did not r e p o r t t h e model used f o r t h e monthly averaged concentration of SOZ.
For illustrative r e a s o n s both DOPT and DOPTEX programs were used to con- s t r u c t the optimal allocation of observation stations and naturally t h e y led t o t h e same (up t o computational a c c u r a c y ) r e s u l t s .
The optimal network consists of seven stations ( t h e model contains six unk- nown paramelers). Usually t h e number of observing stations i s equal to t h e number of unknown p a r a m e t e r s . The seventh point a p p e a r s h e r e due Lo some pecu- l i a r i t i e s in t h e controllable region.
One c a n s e e t h a t t h e v a r i a n c e s of a l l p a r a m e t e r s ( e x c e p t t h e f i r s t one whose v a r i a n c e does not depend upon t h e allocation of stations) are reduced 10-20 times.
Theoretically t h e optimal design assumes t h a t t h e a c c u r a c y of observations at t h e various points i s different. Sometimes t h i s demand i s not r e a l i s t i c in p r a c t i c e but it i s easy t o verify theoretically t h a t t h e design c h a r a c t e r i s t i c s a r e quite s t a b l e under variation of weights (see Fedorov and Uspensky, 1975, p.56). The cal- culations confirm t h i s f a c t f o r our example. For instance, from t h e optimal design ( s e e Printout I), point 1 with small weight ( "0.054) was removed from the design and f o r all o t h e r s t h e weights were chosen equal 1 / 6 (so called s a t u r a t e d design:
number of observation
=
number of unknown parameters). The r a t i o of t h e d e t e r - minants of t h e variance-covariance matrix f o r t h e newly c o n s t r u c t e d design and D-optimal designs w a s found t o b e equal to q , 2 . In terms of variances, t h e discrepancy ( w p m ) i s negligible.Figure 1. Monthly a v e r a g e (January, 1981) of SO2 (in 0.1 ppm) a t 9am f o r Taipei City.
M...
. . ...
YjCI Y......
I I YIR...
I mMS... .em..
...
.an I x !n.... a
.
.a,. In.. .st-n..
... . . . . ..
.nY..
.... ... .
.rVCIYn...n
-.
.nM... M
n
...
M-.I
...
u..
..
.Y(b) Removal of t h e "Least i n j b r m a t t v e " s t a t i o n s . In p r a c t i c e o n e may need to r e d u c e t h e n u m b e r of o b s e r v a t i o n s t a t i o n s avoiding a n y r e l o c a t i o n s . In t h i s case t h e deletion p r o c e d u r e h a s t o b e used; i t c a n b e e x e c u t e d by e i t h e r p r o g r a m . H e r e , DOPTEX i s p r e f e r a b l e from t h e computational point of view.
To b e s p e c i f i c , t h e necessity of removing f o u r s t a t i o n s was assumed. S t a t i o n s 4, 5, 6, 9 (compare Figures 2 a n d 4) w e r e subsequently deleted. The non-normalized d e t e r m i n a n t of t h e v a r i a n c e - c o v a r i a n c e matrix d e c r e a s e d 5.5 times b u t t h e nor- malized m a t r i x i n c r e a s e d 2.7 times. In o t h e r words, t h e d e c r e a s e in t h e non- normalized d e t e r m i n a n t confirms t h a t a n y deletion of d a t a (assuming t h e i r c o r r e c t - n e s s ) d e c r e a s e s t h e final information b u t t h e i n c r e a s e in normalized d e t e r m i n a n t indicates t h a t t h e d e l e t i o n was done in t h e s e n s e of t h e c h o s e n c r i t e r i o n optimally:
t h e "effectiveness" of t h e rest of t h e o b s e r v a t i o n s t a t i o n s significantly i n c r e a s e d ( s e e p r i n t o u t 2).
(c) A d d i t i o n of n e w s t a t i o n s . The "inclusion" algorithm w a s used to find t h e location of t h r e e new s t a t i o n s
-
t h e y a p p e a r e d o t h e boundary of t h e region. The4 "5
non-normalize d e t e r m i n a n t i n c r e a s e d 2.3
-
1 0 times and t h e normalized o n e i n c r e a s e d 5.10 times, i.e., in t h i s case both t h e t o t a l information a n d i t s effective- n e s s i n c r e a s e d ( s e e p r i n t o u t 3 a n d Figure 5).(d) m t i m a l o b s e r v a t i o n n e t w o r k c o n t a i n i n p some s t a t i o n s with fized p o s i - t i o n s . When c r e a t i n g a new o b s e r v a t i o n network, o n e c a n f a c e t h e necessity of including in i t some No ( f o r i n s t a n c e , well equipped) existing stations. If t h e t o t a l number N of s t a t i o n s i s given, t h e n o n e h a s t o c o n s i d e r t h e following design p r o b - lem
€,'=Arg min @[(I-No/N)€
+
(N,/N)t01,€
w h e r e
to
d e s c r i b e s t h e location a n d a c c u r a c y of a n existing s t a t i o n r e q u i r e d t o b e in t h e planned n e t w o r k . The solution of (10) c a n b e computed with t h e h e l p of p r o - g r a m s DOPTEX a n d LINEX.The r e s u l t s of t h e c a l c u l a t i o n s f o r D - c r i t e r i o n a n d model (9) are p r e s e n t e d in P r i n t o u t 4 a n d F i g u r e 6.
E z a m p l e 2. L i n e a r p a r a m e t r i z a t i o n , A - c r i t e r i o n . Theoretically the optimal location of o b s e r v a t i o n a l s t a t i o n s d e p e n d s upon t h e c h o s e n c r i t e r i o n of optimality.
In p r a c t i c e t h e d e p e n d e n c e is usually negligible. To confirm t h i s f a c t , let u s con- s i d e r t h e A - c r i t e r i o n when t h e quality of a location i s c h a r a c t e r i z e d by t h e a v e r -
n
a g e v a r i a n c e of t h e p a r a m e t e r estimators: 4
=
m C D,,=
m -ltr D . Thea =I
r e s u l t s of t h e c a l c u l a t i o n ( p r o g r a m LINOPT) f o r model (9) are p r e s e n t e d in Prin- t o u t 5. The allocation of a l l o b s e r v a t i o n s t a t i o n s coincides (up to computation a c c u r a c y ) with t h a t f o r t h e D-optimal allocation, see Figure 3. T h e m a j o r trace- a b l e d i f f e r e n c e i s in t h e "weights": t h e points which are c l o s e r to t h e o r i g i n h a v e t h e g r e a t e r weights (i.e. t h e a c c u r a c y ( o r number of r e p e t i t i o n s ) of o b s e r v a t i o n s h a s to b e greater f o r t h e " c e n t r a l points").
E z a m p l e 3. N o n l i n e a r p a r a m e t r i z a t i o n , D - c r i t e r i o n .
All p r e v i o u s c o n s i d e r a t i o n s w e r e based o n equation (9) which c a n n o t d e s c r i b e a l l d e t a i l s of t h e p a t t e r n s p r e s e n t e d in Figure 1. Unfortunately Modak a n d Lohani, 1985, did n o t r e p o r t t h e model used to c o n s t r u c t i s o p l e t h s drawn on t h e i r f i g u r e . T h e r e f o r e , w e a p p l i e d t h e following nonlinear r e s p o n s e function t r y i n g t o a p p r o x i - mate t h o s e isopleths:
0.17!x
. . .
x x 0.12 ! x0.08 ! X X . .
. . .
x0.04 !
-0.00 !
. . .
- 1 0 . .-0.04 !
-0.06 !
. . .
x-0.1: !
-0.17 !
. . .
1 I . . . . . .-0.21 !
-0.25 !
. . . .
.l!-0.29 !
-0.3; !
I . . . . .
x-0.36 ! X
-0.42 ! . 2
. . .
x x x-0.46 ! x
-0.50 ! . l x . . . x x x x x x
-0.54 !
-0.5% ! x x x . . .
-0.63 !
-0.67 !
. . .
-0.71 !
-0.35 ! x x
. . .
-0.79 !
-0.23 ! ~8 . . . i
-0.88 ! x
-0.92 ! 1 . 1 . .
-0.96 ! x
-1.00 ! X X X
...
1Xmin
=
-1.000 Xmax = 1.000x - boundary points, O - designs paint
Figure 4: The observation network after deletion of four stations.
Neu p o i n t s appeared on t h e boundary o f the region
1.00 ! x
0.96 !
0.92 ! x . x H
0.S7 !
0.83 ! x x
. . . .
0.75 !
0.25 ! x x
. . .
0.71 !
0.67 ! . . .
0.62 !
. . .
0.56 ! x . . .
0.54 !
0.50 ! x . . O . . . x
0.46 !
0 . 4 2 ! x
. . .
x 0.37 !0.33 !N
. . .
x 0.23 !0.25 !.
. . .
0.21 !
.
0.17!x
. . .
x x0.12 ! x
0.08 ! x x
. . .
0. . .
x 0.04 !x
. . . o . . .
x. 0
. . .
X x x x-0.58 ! x x
I . . .
-0.63 !
-0.67 !
. . .
-0.71 !
-0.75 ! 1 1
. . .
-0.79 !
-0.63 ! x 0
. . .
x-0.88 ! x
-0.92 ! 1 . 1 . .
-0.96 ! x
-1.00 ! x x x
...
Xmin = -1.000 Xmax
=
1.000 x - boundary p o i n t s , 0-
design p o i n t s(they are marked ty 'N')
Figure 5: The observation network with three new stations.
1.00 ! 0 0.96 !
0.92 ! 1 . 1 1
0.87 !
0.83 ! x x . . . .
0.79 !
0.75 ! 1 1
. . .
0.71 !
0.67 !
. . .
0.62 !
0.58 ! I
. . .
0.54 !
0.50 ! x
. .
F . . . x 0.46 !0 . 4 2 ! I . .
. . .
x 0.33 !0.33 !O
. . .
x 0.29 !0.25 !.
. . .
0.21 !
.
0 . 1 7 ! x
. . .
x x0.12 ! x
0.06 ! I X .
. . .
x0.04 !
-0.00 !
. . .
F . .-0.04 !
-0.08 !
. . .
I-0.13 !
-0.17 !
. . .
-0.21 !
-0.25 ! . . .
-0.2Q !
-0.33 ! x
. . .
F . . . x-0.38 ! x
-0.42 !
. . .
x x 0-0.46 ! x
-0.50 ! . F x
. . .
x x x 1 x 1O x x . . . . .
Xain 1 -1.000 Yaax = 1.000
x - boundary points, 0 - neu design points, F - fixed points
Figure 6: D-optimal observation network with five fixed stations.
where aj d e s c r i b e s t h e location of maximums.
The l e a s t s q u a r e method w a s used to f i t t h i s response. The values of concen- t r a t i o n s corresponding to t h e uniform g r i d were t a k e n as "observations" to restore 19'
=
(gll,. . .
,. . . . . .
, 1 9 ~ ~ ) . The corresponding "estimates"( f o r normalized zi : -1 5 zi 5 1 ; i 4 2 ) were found to b e
aT =
(2.83, 4.42, -1.78, 15.3, 1.37, 0.54, -0.50, -2.00). Only t h e f i r s t exponent h a s a bell s h a p e , while t h e second contains t h e negative coefficient I942 corresponding to t h e q u a d r a t i c term;t h i s f a c t tends to confirm t h a t t h e r e a r e more t h a n t w o pollution emission c e n t r e s . One c a n see t h a t in t h i s example, unlike t h e l i n e a r c a s e , w e are concerned with t h e values of t h e p a r a m e t e r s ' estimates. The r e a s o n is t h a t in t h e l i n e a r case t h e variance-covariance matrix d o e s not depend upon estimated p a r a m e t e r s while in t h e nonlinear c a s e ( s e e Fedorov and Uspensky, 1975) this matrix ( o r more accu- r a t e l y i t s asymptotic value) depends upon t h e t r u e values of t h e unknown parame- ters dt :
aa77(2,*) where M(I9.t)
=
Jf (19,zlf ' ( ~ . z ) # ( d l ) , f (I9.z)- aJ .
where N i s t h e number of o b s e r v a t i o n s and
#
is a limit point. Optimal designs for- mally defined by (2) will also depend upon d t , which are naturally unknown a p r i o r i . In t h i s situation t h e following p r o c e d u r e i s recommended:-
a u s e r h a s to choose some p r o b a b l e (reasonable, admissible, e t c . ) values of 19 and define i n t e r v a l s which will almost c e r t a i n l y contain true values of unk- nown p a r a m e t e r s ;-
f o r boundary points of t h e s e intervals, optimal designs have to b e computed with t h e h e l p of one of t h e a b o v e d e s c r i b e d programs;-
if t h e c o r r e s p o n d i n g designs d i f f e r g r e a t l y from e a c h o t h e r , a n "average"design h a s to b e constructed.
Fortunately optimal designs are r a t h e r s t a b l e to t h e variation of p a r a m e t e r s a n d t h e r e f o r e t h e latter p r o c e d u r e c a n b e avoided.
In o u r example, tke v e c t o r
i
w a s t a k e n as a c e n t r a l point and i n t e r v a l s were t a k e n equal to*
0.1-3,. P r i n t o u t 6 and Figure 7 contain information on t h e D-optimal design #f f o r t h e c e n t r a l point. All designs (optimal allocations) practi- cally coincide with
#;
and only in some of them d o one or t w o additional points a p p e a r with small weights. These additions were removed and subsequently t h e corresponding determinants were computed. The r a t i o of determinants f o r t h e modified designs and optimal design fluctuated between 1.02 and 1.09. That i s negligible for p r a c t i c a l needs. T h e r e f o r e ,#f
c a n b e used as a design, defining a n optimal o b s e r v a t i o n network for t h e nonlinear model (11).1.00 ! X 0.Pb !
0.92 ! X . X X
0.27 !
. . .
0.65 ! X X .
0.79 !
0.75 ! ~ J I. . .
0.71 !
0.67 !
. . .
0.62 !
0.58 ! I
. . .
0.54 !
0.50 ! x
. . .
0 0.46 !0.33!x
. . .
x 0.29 !0.25 !.
. . .
0.21 ! .
0.17!O
. . .
O . . . x x 0.12 ! 10.0'2 ! X X . . .
0.04 !
-0.00 !
. . .
-0.04 !
-0.08 !
. . .
0. . .
-0.1: !
-0.17 !
. . .
-0.21 !
-0.25 ! . 0
. . .
-0.29 !
-0.33 !
x . . . . . .
-0 .3b 7 1 I .
-0.42 !
. . .
O . . . -0.46 !-0.50 ! . . I .
. . .
I X-0.54 !
-0.58 ! O X
x . . .
. . . -0.63 !-0.67 ! . . .
-0.71 !
-0.75 ! X I .
. . . .
-0.79 !
-0.83 ! X . . . . I
-0.88 ! x
-0.92 ! X . X . .
-O.?b ! X
-1.00 ! 0 X I !
--- Xmin = -1.000 X a a ~ l- 1.000
x - boundary p o i n t s , 0 - design p o i n t s
F l g u r e 7: D - o p t i m a l observation n e t w o r k f o r nonlinear m o d e l (11).
REFEEZENCES
Atkinson, A.C. & V.V. Fedorov (1987) Optimum Design and Experiments. Encyclo- pedia of S t a t i s t i c a l Sciences. N e w York: Wiley & Sons (to b e published).
Fedorov, V. V. (1972) Theory of Optimum Experiments. New York: Academic P r e s s , pp.292.
Fedorov, V.V. & A.B. Uspensky (1975) Numerical Aspects of Design and Analysis of Experiments. Moscow, Moscow S t a t e University, p.167.
Fedorov, V.V. (1986) Optimal Design of Experiments: Numerical Methods. WP-86- 55. Laxenburg, Austria, International Institute f o r Applied System Analysis Intergovermmental Working Group (1971) P r o c e e d i n g s of Meeting of Intergovern-
mental Working Group on Environmental Monitoring (IGNG, Geneva).
I z m e l , Yu.A. (1974) Global observation system. P r e d i c t i o n and assessment of changes in environmental state. Based o n monitoring. Meteorology: Hydrol- ogy, 7, 3-8 (in Russian).
I z m e l , Yu.A, F.Yu. Rovinsky, M.Yu. Antonovski et al (1984) On t h e s t a t i s t i c a l sub- stantiation of regional atmospheric pollution components in a background area. DAN SSSR 276(2):334-337 (in Russian).
Hanna, S.R., G.A. Briggs & R.P. Hosker (1982) Handbook o n Atmospheric Diffusion.
Technical Information Center, U.S. Department of Energy, p.102.
Mitchell, T.J. (1974) An Algorithm f o r Construction of D-Optimal Experimental Designs. Technomettics, 16203-210.
Mehra, R.K. & D.G. Lainiotis (1976) Systems Identification: Advances and Case Stu- dies. New York: Academic P r e s s , p.593.
Modak, P.M. & B.N. Lohani (1985) Optimization of Ambient Air Quality Monitoring Networks. E n v i r o n m e n t a i Monitoring and Assessment, P a r t 1-11, 5:l-53.
Rovinsky, F.Ya, and G.B. Wiersma (1986) P r o c e d u r e s and Methods of I n t e g r a t e d Global Background Monitoring of Environmental Pollution, WMO. (Forthcom- i w . 1
Silvey, S.D. (1980) Optimal Design. London: Chapman a n d Hall, p.86.
USER'S GUIDE
INSTRUCTIOWS FOR PROGRAM HAP
-
HAPPING OF A CONTROLLABLE REGION R(X)...
---
! SCREEN !
---
I. SPACE DIHEWSIOW
-
? (L) L i s a nurber o f controllable variables 2. Xl(rin), XI(MX)-
?( X l r i n , Xlrax)
Xlrin, Xlmx are the r i a i u l and maximal values o f the f i r s t coordinate 3. GRID FOR X I ? ( ~ 1 ) I n t e r v a l ( Xl(min), Xl(rax)
1
i s divided i n t o NX1 parts, r x defines an i n i t i a l g r i d for XI:
rx = ( X l ( n x )
-
Xl(rin) )/HX1 Hessages 4 - 7 appear i f L-
2.
5. GRID FOB X2 ? (NX2)
XPrin, X2rax
-
minimal and raximal values o f the second coardinate I n t e r v a l ( X2(min), X2(rax)1
i s divided i n t o llX2 parts, r y defines an i n i t i a l g r i d for 12:
r y : ( XP(rar)
-
X2(rin) )/N%zHessage 6 appears for a l l X1
=
x,
belonging t o the grid.6. X I
-
x, BOMB FOR X2 ? Y1 and Y2 are bounds o f the 2-nd(Yl, Y2) coordinate for current value x
of the I - s t coordinate
7. HEM BOUIIDS FOR X2 : I M Y = 1
-
to (6) v i t h the same yes-
1, no-
0 ( I K Y ) value x i f for a given xthe set R
?"
x) i s not come:, R(x)-
{ y : a pair (x,y) belongst o the controllable region )
1
IllEY
=
0-
go t o (6) r i t h ner value x,
x=
x ( m )=
:(old) 4 r xINSTRUCTIOHS FOR PROGRAH DOPT
- ...
FIRST
-
ORDER OPTIHIZATION ALGORITHW FOR D-
CRITERION---
! SCREEN !
---
! COHHEHTS !
1. SPACE DIWENSIOH
-
? (L) L i s a number o f controllable variables 2. COHSTANT FOR CONVERGENCE EPS - a constant for testing con-CRITERIOH
-
? (EPS) vergence o f the algorithm 3. HUHBER OF ESTIHATED PARA- H-
nunber of parameters ( t! r u s tHETERS
-
? (It) correspond t o subroutine RESP,
where a response function i s c a l - culated )
4. HUHBER OF POINTS IH NO
-
number o f supporting points i n INITIAL .DESIGN-
? (HO) an i n i t i a l design5. DESIGN INPUT: FILE
-
1, IDES : I-
i n i t i a l design i s saved MINITOR - 2 (IDES) i n the f i l e 'DES.DAT1;IDES : 2 - i n i t i a l design i s defined on the screen
Hessages 6,7 appear i f IDES = 2 ( i
=
1 ,....,NO )6 . Point i
,
coordinates-
? X(i,k)-
coordinates o f an i - t h point ( X(i,kl, k : 1,L1
i n an i n i t i a l design7. Weight for point i ? P(i1
-
weight of an i - t h point i P ( i ) )6. GRAPHICAL PRESENTATIOH OF INITIAL DESIGH: yes
-
I,no - 0 (ID01
ID0
=
I-
subroutine GRAPH i s executed for i n i t i a l designMessage 9 appears i f i n i t i a l information matrix i s singular.
9. SIKULAR INFORHATIOH HATRIX: IDD : 1
-
go t o ( 4 ) (new i n i t i a l NEW ATTEHPT; yes-
1, no-
0 design i s formed)(1DD) IDD
=
0-
STOP10. SELECTIOH OF 6AIH SEQUEIICE: IALF : 1
-
gain sequence i s constant 1 - a l f a ( s )=
const,
IALF=
2 - gain sequence i s 11s ; 2 - alfa(s) = 11s,
IALF : 3-
steepest descent method3 - alfa(s1 i s steepest i s executed
descent sequence (IALF)
kssage 11 appears i f I A L F
=
I.
11. CONSTAHT FOR 6 A I N SEQUENCE ALFA i s the chosen constant for gain (0.01 - 0.05) (ALFA) sequence ( 0.01 - 0.05
-
recom-mended values
12. NUnflER OF ITERATIOHS ? (NITER)
1;. CONSTAHT FOR HERGIN6 Of SUPPURTING POINTS
-
? (CHER)14. FORYARD LENGTH Of EXCURSION (IIFOR)
15. BACKYARD LENGTH OF EXCURSION ( n a A c )
16. I N I T I A L PROCEDURE:
forward
-
1, backward-
2 (IPRO)17. STEPYIZE IHFURHATIOH :
yes - 1, no
-
0 ( I I N F )NITER - maximal nurber o f i t e r a t i o n s
CHER i s an i n t e r n a l constant ( see section 3
1
NFOR
-
number of steps for forward procedureNBAC - number o f steps f o r backward procedure
The algorithm s t a r t s with:
- forward procedure i f IPRO = 1.
-
backward procedure i f IPRO-
2.I I W F
-
1-
intermediate information i s saved i n the f i l e '0UT.DAT' and shown on the monitor (current design, value o f the determinant etc)Hessage 18 appears i f L
=
2.
18. GRAPHICAL PRESENTATION OF I 6 R : 1
-
subroutine GRAPH i s DfSI6II: yes-
I, no - 0 executed for f i n a l design( I 6 R )
19. SCALIHG OF DESIGN: yes - 1, ISC : 1
-
scaling o f f i n a l designno
-
0 ( I S C ) i s c a r r i e d outHessages 20
-
22 appear i f I S C : 1.20. Xl(min), Xl(max)
-
? (Xlmin, Xlmax)Hessage 21 appears i f L
=
2.
22. GRAPH I N REAL SCALE:
yes
-
1, no-
0 (IGRS)Xlmin, Xlmax - minimal and maximal values o f the 1-st coordinate
XPmin, X2max
-
minimal and maximal values o f the 2-nd coordinateI G R S
=
I-
subroutine GRAPH i s executed for f i n a l design i n r e a l scale.INSTRUCTIOHS FOR PRLiGRAH DOPTEX
-
O P T I H I Z A T I O H ALGORITHH Of THE EXCHANGE TYPE FOR D
-
CRITERION...
--- --
! SCREEN !
---
---
! COHHENTS !
---
Hessages 1
-
9 coincide with those for program DOPT.
10. HUHBER OF F I X E D POINTS I N 'The first H F I X points in initial &sign I H I T I A L DESIGN ( K I X ) are fixed
11. CHOICE OF THE ALGORI'THH: I A L G
=
1-
deleting (ONLY !) is1
-
DELETING PROCEDURE, executed;2
-
INCLUDIHG PROCEDURE, I A L G=
2-
including (ONLY !) is3
-
STANDARD EXCHANGE PRO- e~ecuted;CEDURE ( I A L 6 ) IM6
=
3 - exchange procedure with including and deleting is executedHessage 1 2 appears if I A L 6 1 1.
12. NUHBER OF POINTS FOR NBAC
-
number of steps for deleting DELETING (NBAC) procedure ( M A C mst be less thanHO
-
H !!! )Message 13 appears if I A L G : 2.
13. HUHBER OF POXHTS FOR NFOR
-
nurber of steps for includingI N C L U D I K (NFOR) procedure
Hessages 1 4
-
19 appear if I A L G = 3.1 4 . SELECTION OF STEPSIZE SEQUEHCE:
1
-
a l f a = const, 2 - a l f a=
11s ( I A L F )Message 1 5 appear if I A L F = 1.
15. COHSTANT FOR STEPSIZE
-
? (ALFA)I A L F
=
1-
stepsize is constant;I K F
--
2-
stepsize is of the form 11s,
s=
HO+l, NOt2,...
ALFA is a constant stepsize for the algori thr
16. NUHEER OF ITERATIONS
-
? HITER-
maximal number of i t e r a t i o n s (MITER)17. CONSTANT FOR HERGING OF CHER. i s an i n t e r n a l constant o f the SUPPORTING POIMTS (CKR) a l ~ o r i t h m (see scctinn 3 )
18. LENGTH OF EXCURSION - ? (forward and backward)
(NFOR)
NFOR - n u d e r of 'steps for forward and backward procedures
( Attention :
HITER = 2*NFORtK. K
-
i n t e j e r ! !!) Messages 19 - 25 coincide with messages 16 - 22 forprogram DUPT
.
Instructions f o r programs LINOPT and L I K X are almost the same as for programs DOPT and DOPTEX respectively. There appears one additional ressage ( a f t e r message 3 1
.
3*. UTILITY HATRIX { u t ( i , j ) 1 i s a symmetric u t i l i t y
(upper triangular p a r t ) .
.
matrix,( ut(lI1), ut{1,2),
...
ut(1,m) J : 1,...
, m i i=
I , . - . , m . ut(2,2),... ,
ut(2,m)...
ut(m,m)
INSTRUCTIOMS FUR SUBROUTINE GRAPH
-
GRAPHICAL PRESENTATIOH OF THE DESIGN
! SCREEN !
---
I. Hurber o f divisions for X I ? (HX)
2. Number o f divisions for X2 ? (HY1
The graph has HX positions for the f i r s t coordinate
and
HY positions f o r the second coordinate
-
27-
Printout 1: Completely new network, D-criterion.
Exmple for prograr DCPT: Oljtpilt f i l e 'QUT.DAT'
SPACE DiHENSIuN - ? 2
CI)HSTA)IT FOR CONVERGENCE CRITERION - ? 0.01000
NUHGER OF ESTIHATED FARAHETERS - ? 6
NUHRER OF POINTS I N INITIAL DESIGN - 7 11
INITIAL DESIGN p o i n t v e i g h t 1. 0.091 2. 0.05'1 3. 0.091
4 . 0.091
5. 0.0?1 6. 0.091 7. 0.091 8. 0.091 9. 0.0?1 10. 0.091 11. 0.0?1
t t S S S S t S S
coordinates -0.5789 -0.5000 -0.5369 -0.4167 -0.4737 0.5000 -0.3663 -0.5000 -0.3684 -0.166;
-0.2632 0 . 0 l j 3 -0.1579 -0.1667 -0.0526 -0.8333 0.1579 -0.3333 0.3624 0.
0.4737 -0.2500
INITIAL INFORHATION HATRIX 1
.
000-0.167 0.150 0.150 -0.041 0.035 -0.235 0.038 -0.032 0.189
0.169 -0.042 0.023 -0.076 0.065 0.056 -0.032 0.012 -0.032 0.011 0.023
DETERHIWANT OF INITIAL IHFORHATIOH HATRIX 8.08280e-09
INITIAL COVARIANCE HATRIX 6.135
5.457 30.835
-21.406 -7.066 123.469 13.226 59.169 - 3 . 8 3 1 85.518
5.754 36.S72 -23.188 74.240 96.OE7
29.515 89.311 -38.266 167.743 144.597 401.929
Number o f c o n t r o l l a b l e variables :
L . 2
Constant f o r convgrgence test: eps = 0.01
Hurber o f parameters : n = d
Humber o f p o i n t s i n i n i t i a l design : NO : 11
Supporting p o i n t s o f the i n i t i a l design and t h e i r weights
I n i t i a l information m a t r i ~ : since the matrix i s s y r a e t r i c , only i t s l o u t r i a n g u l a r p a r t
i s shown
The value o f the o p t i m a l i t y c r i t e r i o n : i n t h i s example
i t i s the deterainant o f the i n f o r n a t i o n matrix
I n i t i a l covariance matrix:
for example, d(21 = 30.235 i s the variance o f the i n i t i a l estiaator f o r the 2-nd parameter
