Velten K. Mathematical Modeling and Simulation: Introduction for Scientists and Engineers

Подождите немного. Документ загружается.

2.6 Design of Experiments 105

of these locations. But then, the location itself may affect the hardness reading.

Referring to the metal test coupon ‘‘A’’ in Figure 2.15, for example, the hardness

measurements using tips 2 and 3 may be somewhat inﬂuenced by the fact that their

measurement positions are more close to the ends of the metal coupon compared

to the measurement positions of tips 1 and 4.

If the hardness measurement really depends on the number of measurements

that have been previously performed on the same metal test coupon, or if the

measurement depends on the measurement location on the metal test coupon, a

designsuchastheoneshowninFigure2.15–whereoneofthetipsisusedmore

than once in one particular position of the measurement order (such as tip 4) – is

obviously an unsuitable design. Based on such a design, differences between the

hardness readings produced by the tips may be observed which are caused by

the particular position of the measurement order where the tips are used (rather

than by the tips itself). Obviously, a better experimental design should randomly

distribute the factor levels on each of the metal test coupons in a way such that any

factor level appears only once in one particular position on the metal test coupons.

A design of this kind is known as a Latin square design.

In R, a Latin square design for the above example can be computed using the

following code (see

LSD.r in the book software):

1: library(agricolae)

2: levels=c("Tip 1", "Tip 2", "Tip 3", "Tip 4")

3: out=design.lsd(levels,number=1)

5: print(out)

(2.80)

This code may yield the following result (may: see the discussion of code (2.6.1)

in Section 2.6.1):

plots row col levels plots row col levels

1 111Tip1 9 931Tip4

2 212Tip2 10 1032Tip1

3 313Tip3 11 1133Tip2

4 414Tip4 12 1234Tip3

5 521Tip3 13 1341Tip2

6 622Tip4 14 1442Tip3

7 723Tip1 15 1543Tip4

8 824Tip2 16 1644Tip1

Figure 2.16 visualizes this result similar to Figure 2.15 above. As can be seen,

every tip occurs only once in each single row and column of this experimental

design as required.

The design discussed so far can be generalized in various ways. While the Latin

square can be used to treat situations with two sources of extraneous variability

(in the above example, the metal test coupon and the position of a tip in the test

sequence), the Graeco–Latin square design can treat a similar situation with three

106 2 Phenomenological Models

Fig. 2.16 Latin square design computed using design.lsd

in R (hardness testing example): metal test coupons A–D,

and numbers indicating the randomly chosen sequences of

the tips 1–4.

sources of variability. In R’s agricolae package, Graeco–Latin square designs

can be used based on the

design.graeco command in a similar way as the

corresponding commands that were discussed above. In some cases, a RCBD may

be too demanding in terms of time and resources, that is, we may be unable to

test every factor level in each block of the design. Then, a randomized balanced

incomplete block design can be used, which can be obtained using the

design.bib

function of R’s agricolae package.

2.6.4

Factorial Designs

All experimental designs considered so far involved one factor only. If an experiment

involves two factors or more, factorial designs or response surface designs are applied

[73]. Factorial designs are preferentially used in situations where each factor is

varied on two levels only, a lower level which is typically designated as ‘‘−’’,

and an upper level designated as ‘‘+’’. Designs of this type are called two-level

factorial designs. Response surface methods, on the other hand, focus on situations

where the factors are varied on more than two levels. They provide procedures

that can be used to decide about an optimal choice of the factor levels, based

on an approximation of the response of the system (y) depending on the factors

, x

, ...) e.g. using polynomials. We will conﬁne ourselves here to the basics of

factorial designs (see [72, 73] for more on response surface methods).

A factorial design involving n ∈ N factors on m ∈ N levels is usually denoted as a

‘‘m

design’’. Hybrid factorial designs such as a 2 × 3

design are also used, which

involves one factor that is varied on two levels and two factors that are varied on three

levels in this case. As was already mentioned, the most frequently used factorial de-

signs are 2

designs, that is, designs involving n factors that are varied on two levels.

As a simple example, consider a chemical reactor which is used to produce some

product, and which is affected by two factors: the amount of a catalyst that is used

and the temperature. Let us denote the catalyst with A and the temperature with B,

and let us assume that an experiment is performed where both factors are varied

at a low and a high level, respectively, which we denote as ‘‘−’’ and ‘‘+’’. Table 2.2

2.6 Design of Experiments 107

Table 2.2

Example of a 2

factorial design.

No. A B Yield

1 −− 50

2 +− 54

3 −+ 64

4 ++ 90

shows a possible result of such an experiment (the yield of the product is given

in appropriate units which we do not need to discuss here). In the terminology

introduced above, this is a 2

design as it involves two factors being varied on two

levels. Experimental designs that involve all possible combinations of the factor

levels, such as the design in Table 2.2, are also called full factorial designs.

Full factorial designs such as the design in Table 2.2 can be generated very easily

using R’s

expand.grid command as follows:

1: levels=c("-", "”)+

2: expand.grid(A=levels,B=levels)

(2.81)

which yields the design that was used in Table 2.2:

AB AB

1-- 3-+

2+- 4++

For the same reasons that were explained in Section 2.6.2 above, the exper-

imenter may decide to use a randomized block design in a factorial experiment.

In R,suchadesigncanberealizede.g.usingthe

design.ab command of R’s

agricolae package. The package contains an example experiment that involves

the ‘‘perricholi’’, ‘‘canchan’’, and ‘‘tomasa’’ potato varieties which are cultivated

using three nitrogen levels. A full factorial 3

design involving ﬁve replications of

each experiment is used, and the replications are organized in randomized blocks.

To compute such a design, the following code can be used in R (see

FacBlock.r

in the book software):

1: library(agricolae)

2: variety=c("perricholi", "canchan", "tomasa")

3: nitrogen=c(40,80,120)

4: out=design.ab(variety, nitrogen, 5, number=1)

5: print(out)

(2.82)

This code works very similar to the codes discussed in the previous section. Since

we use a 3

design here which involves ﬁve replications, you can easily anticipate

108 2 Phenomenological Models

that this code will result in a total of 3

· 5 = 45 experimental runs. The ﬁrst part of

theoutputlookslikethis:

plots block variety nitrogen plots block variety nitrogen

1 1 1 tomasa 80 10 10 2 canchan 80

2 2 1 perricholi 120 11 11 2 canchan 40

3 3 1 canchan 120 12 12 2 perricholi 120

4 4 1 perricholi 80 13 13 2 tomasa 80

5 5 1 canchan 80 14 14 2 perricholi 40

6 6 1 perricholi 40 15 15 2 canchan 120

7 7 1 tomasa 40 16 16 2 perricholi 80

8 8 1 canchan 40 17 17 2 tomasa 40

9 9 1 tomasa 120 18 18 2 tomasa 120 ...

As can be seen, the code produces a complete 3

design involving all possible

combinations of the potato varieties and nitrogen levels (randomly ordered) in each

of the blocks.

Of course, time and resources restrictions may prevent us from performing full

factorial experiments involving 45 different experimental runs, and similar to the

‘‘randomized balanced incomplete block design’’ discussed in Section 2.6.3 above,

the experimenter will be interested in clever ways to reduce the overall number

of experimental runs in a way that does not affect the signiﬁcance of his results.

Factorial designs that do not consider all possible combinations of the factor levels

are called fractional factorial designs. Simple fractional factorial designs are usually

denoted as ‘‘m

n−k

designs’’, where m (the number of factor levels) and n (the

number of factors) have the same meaning as before, and k expresses the fact

that the total number of experimental runs has been reduced by a factor of 1/2

Example: While a full factorial 2

design would require 2

= 32 runs, a fractional

5−2

design needs only 2

= 8 runs. To generate such fractional factorial designs,

one can use R’s

ffDesMatrix command (a part of the BHH2 contributed package).

See [72, 73] for more theoretical background on factorial and fractional factorial

designs.

2.6.5

Optimal Sample Size

An important issue in experimental design is the proper selection of sample size,

that is, of the number of repetitions of a particular experiment that are necessary to

achieve the desired results. This is very well supported by a number of R functions

such as

power.t.test, power.anova.test, power.prop.test,etc.,whichare

apartof R’s standard distribution. We will conﬁne ourselves here to a simple

example that demonstrates the way in which these functions can be used. Suppose

the yield of crop variety A is suspected to be higher than the yield of crop variety

B, and suppose you want to show this using a one-sided t test with α = 0.1and

2.7 Other Phenomenological Modeling Approaches 109

β = 0.9 (see Section 2.1.3 for details on the testing procedure). Then, the following

R command can be used to compute an optimal sample size:

power.t.test(sig.level=0.1,power=0.9,delta=2,sd=1

,alternative="one.sided")

This yields an optimal sample size of n ≈ 3.9inR, which means that n = 4

replications should be used for each of the crop varieties. The

delta argument is

the ‘‘true difference in means’’, that is, μ

− μ

(see Section 2.1.3). Of course, you do

not know this true difference in means a priori in a practical situation, so it should

be set to a difference that you want to be detected by the experiment. For example,

a difference in the yields of the varieties A and B may be practically important only

if it is larger than 2 kg m

−2

, so in this situation you would set delta=2 as above.

Of course, it is more difﬁcult to get a signiﬁcant test result if

delta is small, which

is reﬂected by the fact that n increases as

delta decreases (for example, n = 14 if

delta=1 is used in the above command). Hence delta should be chosen as large

as possible such that it still satisﬁes the above requirement. The

sd argument of

power.t.test is the standard deviation of the random variables that generate the

data, which is assumed to be constant here across the crop varieties. Again, this is

not known a priori, and in a practical situation you would set

sd according to the

experience made in similar prior experiments (if available), or you would have to

guess an appropriate order of magnitude based on your knowledge of the system

and the measurement procedure.

2.7

Other Phenomenological Modeling Approaches

You should note that there is a great number of phenomenological modeling ap-

proaches beyond the ones introduced above. Only a few of these topics can be brieﬂy

addressed in the following sections: soft computing approaches in Section 2.7.1,

discrete event simulation in Section 2.7.2, and signal processing in Section 2.7.3.

2.7.1

Soft Computing

Soft computing is used as a label for relatively new computational techniques such

as artiﬁcialneuralnetworks(ANNs), fuzzy logic, evolutionary algorithms,butalsofor

recent developments in ﬁelds such as rough sets and probabilistic networks [74, 75]. As

it is explained in [75], a common feature of these techniques is that, unlike conven-

tional algorithms, they are tolerant of imprecision, uncertainty, and partial truth.

Artiﬁcial neural networks have already been introduced in Section 2.5. The

above discussion was restricted to neural networks used as a tool for nonlinear

regression. As mentioned there, there is a great number of applications in various

110 2 Phenomenological Models

other ﬁelds such as time series prediction, classiﬁcation and pattern recognition,

or data processing. An important feature of neural networks is their ability to

‘‘learn’’ from data. We have seen in Section 2.5 that a neural network is able to

‘‘learn’’ the nonlinear shape of a function from a dataset, that is, there is no need

to specify this nonlinear shape as a mathematical function as it is required in the

classical nonlinear regression approach. This kind of adaptivity, that is, the ability

of a model to adapt to a changing problem environment, is a characteristic feature

of soft computing approaches in general [75].

As explained above, artiﬁcial neural networks were originally inspired by an

analogy with biological neural networks (Note 2.5.2). In a similar way, evolutionary

algorithms encompass a class of stochastic optimization algorithms that were

originally inspired by an analogy with the biological ideas of genetic inheritance

and the Darwinian law of the ‘‘survival of the ﬁttest’’. In genetic algorithms –the

most widely used type of evolutionary algorithms – individuals are represented

as arrays of binary digits that can take on the values 0 or 1. Basically, these

arrays can be thought of as representing the genes of the individuals. After a

random initial population has been generated, an iterative process starts where

new generations of the population are generated from the previous population

by applying a certain number of stochastic operators to the previous population,

which basically can be thought of as reﬂecting the Darwinian law of the ‘‘survival

of the ﬁttest’’. Similar to neural networks, this bio-inspired approach turned

out to be extremely fruitful in the applications. Evolutionary algorithms have

been applied in bio–informatics, phylogenetics, computer science, engineering,

economics, chemistry, manufacturing, mathematics, physics, and other ﬁelds. See

the examples in [76–78], and [75] for a detailed case study involving a ﬁnancial

application (portfolio optimization). Note that evolutionary algorithms are a part

of the larger ﬁeld of evolutionary computation which includes other techniques

such as swarm intelligence, which describe the collective behavior of decentralized,

self-organized systems [79]. Evolutionary computation itself is usually classiﬁed as

asubﬁeldof artiﬁcal intelligence, a discipline of computer science.

As regards software for soft computing applications, we have already used

R’s

nnet package in Section 2.5 above to do neural network-based nonlinear

regression. The same package can also be used to solve classiﬁcation problems,

see [45]. R’s contributed package

genalg (R Based Genetic Algorithm) can be used

to implement genetic algorithms. While

nnet comesasastandardpartoftheR

distribution,

genalg can be obtained from R’s internet site (www.r-project.org). R

may also serve as a platform for the implementation of fuzzy models (there are a

number of fuzzy-based contributed packages, see the list on www.r-project.org).

2.7.1.1 Fuzzy Model of a Washing Machine

We end this section on soft computing with an example of a fuzzy model. Based

on the ideas developed in [80], fuzzy models use logical variables that can take on

any value between 0 and 1. In this sense, these models allow for ‘‘uncertainty’’,

as opposed to the usual concept where logical variables can take on the values 0

and 1 only. This is of interest in many technological applications where a system

2.7 Other Phenomenological Modeling Approaches 111

needs to be controlled in a smooth way, that is, not based on conventional ‘‘on/off’’

switches, but rather based on a kind of control that allows a smooth transition

between the various states of a system.

An example of a washing machine controlled by a fuzzy model is discussed in

[81]. In that example, the amount of detergent that is used by the machine must

be determined based on the dirtiness of the load (as measured by the opacity of

the washing water using an optical sensor system), and based on the weight of the

laundry (as measured by a pressure sensor system). Using practical experiences, a

set of control rules can be established which determines the amount of detergent

that should be used as the dirtiness varies between the classes of a so-called fuzzy

subset:

Almost_Clean, Dirty, Soiled, Filthy

and as the weight of the laundry varies between the classes of another fuzzy subset:

Very_Light, Light, Heavy, Very_Heavy

Now practical experience may tell us that the amount of detergent should be

increased by a certain amount in a situation where the weight of the laundry is

in the class

Light while the dirtiness increases from Soiled to Filthy. Then, a

fuzzy model will deﬁne a smooth transition of the amount of detergent that is used

as the dirtiness changes its class from

Soiled to Filthy. Basically, it will allow a

classiﬁcation of dirtiness partially between the classes

Soiled and Filthy,thatis,

it will allow for uncertainty as explained above, and the amount of detergent that

is used will reﬂect this uncertainty in the sense that it will be an average amount

between the amounts that are deﬁned in the rules for the classes

Soiled and

Filthy. Control strategies of this kind often produce better results compared to

classical ‘‘on/off’’ strategies. See [75, 81] for a number of other examples such as

vacuum cleaners, antilock brakes, and so on.

Fuzzy models can be described as an approach that leaves the realms of

conventional mathematics to some extent. Beyond fuzzy models, there is a great

variety of other modeling approaches that are located somewhere in a transition

zone between quantitative (i.e. mathematical model based) and qualitative modeling

approaches. See [82] for examples of qualitative approaches used in the social

sciences (such as narrative analysis, action research, critical ethnography, etc.),

and [83, 84] for a case study approach that uses both quantitative and qualitative

methods.

2.7.2

Discrete Event Simulation

In this book, you have already seen (and will continue to see) that there is a great

number of different approaches in mathematical modeling and simulation. With

this in mind, you may be surprised to know that books such as the ‘‘Handbook

112 2 Phenomenological Models

of Simulation’’ by J. Banks (ed.) [6] or ‘‘Modern Simulation and Modeling’’ by

Rubinstein and Melamed [85] both are devoted to one of these approaches only:

discrete event simulation. Two things can be derived from this: (i) These discrete

event simulation people have a very substantial self-esteem and self-conﬁdence,

and they know that what they do is an important part of the overall effort in

modeling and simulation. (ii) You should know what they are doing!

Discrete event simulation can be deﬁned as follows [10]:

Deﬁnition 2.7.1 (Discrete event simulation) Discrete event simulation con-

cerns the modeling of a system as it evolves over time by a representation in

which the state variables changes instantaneously at a (countable number of)

separate points in time. These points in time are the ones at which an event

occurs, where event is deﬁned as an instantaneous occurrence that may change

the state of a system.

To understand the way in which discrete event simulation is used in practice,

consider the following example which we cite here (almost unchanged) from [86]:

Suppose that a single server (such as a clerk, a machine, or a computer) services

randomly arriving customers (people, parts, or jobs). The order of service is ﬁrst

in, ﬁrst out. The time between successive arrivals has the stationary distribution F,

and the service time distribution is G.Attime0,thesystemisempty.

We may simulate this system as follows. First, generate a random number A

using the distribution F,andthenarandomnumberB

using G. See Section 2.1.2.3

for details on distributions, the discussion of

RNumbers.r in Sections 2.1.2.3 and

2.1.2.5 for an example of random number generation using R, and [86] for the

theoretical background of random number generation. Customer 1 enters the

system at time A

(which is an event in the sense of the deﬁnition above) and leaves

at time A

+ B

(another event). Next, generate two more random numbers A

and

. Customer 2 arrives at A

+ A

.IfA

+ A

≥ A

+ B

, he starts service right

away and ﬁnishes at A

+ B

. And so on. A simulation of this kind can be

used to estimate the number of customers that are served up to a given time, their

average waiting time, and so on.

In a similar way, discrete event simulations can be used to simulate a great

number of systems in various ﬁelds. They are used for example, to optimize

and increase the efﬁciency of systems in manufacturing and material handling,

logistics and transportation, and healthcare, see [6, 10]. Discrete event simulations

are supported on a number of commercial and open source software platforms

such as:

•

simcol, open source, contributed package of R,see

www.r-project.org and [87]

•

OpenModelica,opensource,see

www.ida.liu.se/labs/pelab/modelica/OpenModelica.html and [7]

•

jemula, open source, see jemula.origo.ethz.ch/

•

eM-Plant, commercial, see www.ugsplm.com

2.7 Other Phenomenological Modeling Approaches 113

•

SIMPROCESS, commercial, see www.caci.com

•

simul8, commercial, see www.simul8.com/

See Section 4.11.2 for an application of R’s

simcol package.

2.7.3

Signal Processing

As it will become clear in the following chapters on differential equations, many

mathematical models involve rates of changes of quantities of interest (see Note

3.1.1). Hence, it is a natural task to compute rates of changes from experimental

data. As an example, let us reconsider the dataset

spring.csv that was already

analyzed in Section 1.5.6 above. Table 2.3 shows the x and y data of this dataset,

and let us assume that we are interested in the rate of change of y with respect to

x, that is, in the derivative y



(x). Now the question is how y



(x) can be derived from

discrete data of the general form (x

, y

), (x

, y

), ...,(x

, y

). In a naive approach,

we could use the approximation



) ≈

y

x

− y

i−1

− x

i−1

(2.83)

for i = 2, ..., n, which is based on the deﬁnition of the derivative as



(x) = lim

h→0

y(x + h) − y(x )

(2.84)

Using this approximation in spring.csv yields the values labeled as ‘‘y/x’’

in Table 2.3. As can be seen, these values scatter substantially in an interval

between 0.1 and 0.6, that is, it seems that there are substantial changes of y



(x)asx

is increased. Looking at the plot of the data in Figure 1.9 (Section 1.5.6), however,

it seems more likely here that the scattering of the y/x-data computed from the

naive approach is caused by measurement errors only, which make the data swing

a little bit around the regression line in Figure 1.9b that represents the ‘‘true’’

relation between x and y. Thus, we see that the naive approach is very sensitive

to measurement errors, which makes this approach unusable for what is called

the numerical differentiation of data, that is, for the determination of approximate

derivatives from discrete datasets.

A better approach is already suggested by our analysis of

spring.csv in Section

1.5.6 above, where we used the regression line in Figure 1.9b as a mathematical

Table 2.3 Naively computed rates of change in the dataset spring.csv.

x 10 20 30 40 50

y 3 5 11 12 16

y/x 0.2 0.6 0.1 0.4

114 2 Phenomenological Models

model of the data, which was expressed as

y(x) = 0.33x − 0.5 (2.85)

As explained there, this equation expresses the general (linear) tendency of

the data, but it ‘‘damps out’’ the oscillations in the data that are e.g. induced by

measurement errors. Equation 2.85 tells us that the rate of change of y with respect

to x is constant:



(x) = 0.33 (2.86)

This contradicts the conclusion that was drawn above from the naive approach,

but a constant rate of change of y with respect to x obviously is in much better

coincidence with the data that we see in Figure 1.9. The idea of differentiating

regression functions is used quite generally as one of the standard procedures for

the numerical differentiation of data. In the general case, the data may of course

be nonlinear, which means that one will have to use, for example, polynomials as

regression functions instead of the regression line that was used above. To obtain

an approximation of y(x

), one will typically ﬁt a polynomial to a few datapoints

around x

only. This procedure can be viewed as an implementation of certain low

pass ﬁlters, termed variously as Savitzky–Golay smoothing ﬁlters, least squares ﬁlters,

or DISPO (digital smoothing polynomial) ﬁlters [88–90]. See [91] for an example

application of this method to an analysis of the resin transfer molding (RTM)

process, a process that is used in the manufacturing of ﬁber-reinforced composite

materials.

There are many other approaches that can be used to ‘‘damp out’’ measurement

error-induced oscillations of data, such as the moving average approach where a

particular measurement value y

basically is replaced by the average of a certain

number of neighboring measurement values, see [92]. These methods are a

part of the large ﬁeld of signal processing, which is concerned with the analysis,

interpretation, and manipulation of signals, preferentially signals in the form of

sounds, images, biological signals such as the electrocardiogram (ECG), radar

signals, and so on, but its methods are also used for the analysis of general

experimental datasets [93]. Anyone who has used a MP3 player knows about the

power of modern signal processing methods. Signal processing covers approaches

such as the famous Fourier analysis, which allows a decomposition of a function

in terms of sinusoidal functions, and which is used e.g. to remove unwanted

frequencies or artifacts from audio or video recordings, or the recently developed

wavelet transforms, which can be roughly described as a further development of

the general idea of the Fourier transform and which have a similarly broad range

of application (they are particularly well known for their use in data compression

applications).

The open source software R provides a number of signal processing–related

packages. Examples are the

decompose package that can be used for the decompo-

sition of a time series dataset into seasonal, trend, and irregular components using