Scientific and Engineering Applications Using Matlab
Scientific and Engineering Applications Using Matlab
Scientific and Engineering Applications Using Matlab
ENGINEERING
APPLICATIONS
USING MATLAB
Edited by Emilson Pereira Leite
Scientific and Engineering Applications Using MATLAB
Edited by Emilson Pereira Leite
Published by InTech
Janeza Trdine 9, 51000 Rijeka, Croatia
Statements and opinions expressed in the chapters are these of the individual contributors
and not necessarily those of the editors or publisher. No responsibility is accepted
for the accuracy of information contained in the published articles. The publisher
assumes no responsibility for any damage or injury to persons or property arising out
of the use of any materials, instructions, methods or ideas contained in the book.
Scientific and Engineering Applications Using MATLAB, Edited by Emilson Pereira Leite
p. cm.
ISBN 978-953-307-659-1
free online editions of InTech
Books and Journals can be found at
www.intechopen.com
Contents
Preface IX
The purpose of this book is to present 10 scientific and engineering works whose nu-
merical and graphical analysis were all constructed using the power of Matlab® tools.
The first five chapters of this book show applications in seismology, meteorology and
natural environment. Chapters 6 and 7 focus on modeling and simulation of Water
Distribution Networks. Simulation was also applied to study wide area protection for
interconnected power grids (Chapter 8) and performance of conical antennas (Chapter
9). The last chapter deals with depth positioning of underwater robot vehicles. I al-
ways tell my students that by using interactive software like Matlab®, they can focus
their efforts on learning and applying difficult concepts rather than in details of pro-
gramming. Therefore, this book is a collection of interesting examples of where this
computational package can be applied. As the Editor of this book, I would like to
thank InTech - Open Access Publisher for all the support during the publishing pro-
cess and authors for their efforts in order to produce high quality works.
1. Introduction
This chapter presents the implementation of seismological and earthquake engineering
principles and the development of innovative computer code using Matlab platform. Since
earthquakes remain the greatest natural disaster for modern society, causing loss of life and
millions of damages to the urban environment, the efforts of earth scientists and engineers
lie in providing the tools for quick and adequate assessment of potential damage.
The foundation of seismic hazard analysis is based on the accurate scientific estimation of
anticipated ground motion at a site following the occurrence of a strong earthquake. The
seismic parameters involved in this estimation are the magnitude of the earthquake, the
distance between the epicenter and the site in question, the description of the site’s
geological formations and additional characteristics of the earthquake’s rupture style.
Nowadays, knowledge about the level of the anticipated shaking near cities or villages is
directly linked with past observations. In scientific practice the above statement describes
the development of elaborate empirical mathematical models, which are based on the
available seismological data. Seismological data collection and analysis is a demanding
time-consuming task, related with Digital Signal Processing and Data Archiving. The goal of
this chapter is two-fold; firstly to provide an insight of the necessary Digital Signal
Processing steps, easily performed through Matlab, leading to the derivation of earthquake
engineering parameters and secondly testing traditional regression analysis and
optimization in order to develop empirical equations modeling the aforementioned
parameters.
2. Data processing
According to modern data acquisition practice once an earthquake, exceeding a specific
threshold occurs, ground motion time-series recorded by digital seismometers or
accelerometers, usually at a sampling frequency equal to 200 samples-per-second, are
transmitted to the data analysis center. No matter the progress of modern technology,
scientists in earthquake prone countries cannot simply ignore the older analog acceleration
time series, which in many cases can be irreplaceable. Matlab through resample1 function
provides the necessary tool to produce an equally sampled time series at a given sampling
interval. After this processing stage, the main objective is to remove the undesirable long
period and high frequency noise, which can be attributed to various sources like the
2 Scientific and Engineering Applications Using MATLAB
presents a comparison between causal and phase–preserving filtering during strong motion
processing. Figure 3 presents the uncorrected acceleration and filtered acceleration time
series using five different pairs of cut-off frequencies combined in a phase preserving pass-
band Butterworth filter‘s implementation.
Fig. 2. Comparison between phase-preserving (in red) and causal (in blue) implementation
of an IIR filter.
Another important aspect is the integration of a sinusoidal signal, such as the acceleration
time series in this case study (Figure 1). This computation in the frequency domain
corresponds to the convolution of Fourier amplitude spectrum with the frequency response
of the perfect integration operator, equal to 1/iω, whereas for differentiation the frequency
response of the operator is just the inverse of the perfect integrator, simply iω (Karl, 1989).
After convolution the user can easily select the real part of an array of complex numbers,
such as the Fourier amplitude spectrum, by using the real function.
Immediately after filtering the acceleration time series the inspection of velocity and
displacement time series (Figure 4), calculated after single and double integration
respectively, is required in order to determine whether the high-frequency and long period
noise has been removed adequately from the records.
In more elaborate mathematical calculations, related with the response of a single degree of
freedom (SDOF) harmonic oscillator u ɺɺ of a specific damping level ζ subjected to an
acceleration time series x (Equation 1a), the computational effort required is greater. In order
to calculate the response spectral acceleration SA at a given period ω, the user should define
the maximum of the oscillator time series (Equation 1b) for this specific period.
4 Scientific and Engineering Applications Using MATLAB
xɺɺ + 2 ζωxɺ + ω2 x = u
ɺɺ(t )
(1a)
S A ( ωn , ζ ) = max xɺɺ ( t ) + u
ɺɺg ( t ) (1b)
t
The calculation of the response spectral acceleration of the damped SDOF harmonic
oscillator over a range of periods and various damping levels provides the response
acceleration spectrum (Figure 5), which describes the shaking of typical structures during an
earthquake. Figure 5 is the output of Proschema software using the Plot Pseudo Spectral
Acceleration option.
Fig. 3. Uncorrected (in blue) and filtered acceleration (in red) time series.
Ground Motion Estimation During Strong Seismic Events Using Matlab 5
Fig. 5. Response acceleration spectrum for various damping levels (0%, 2%, 5%, 10% and
20% of the critical damping) over two hundred period estimators ranging between 0.01 s
and 10.00 s.
2
log 10 ( PGA ) = a+ bM+ cM 2 + ( d+ eM ) log 10 R 2 + ( H- h ) + e + f (2)
In Equation (2) the seismic magnitude M, distance R (km) and depth H (km) are considered
to be the independent variables of the model (Equation 2). Random variables were
introduced in Equation (2) for modeling soil site conditions (e) and style of faulting (f).
From this point on the scientific effort focuses in solving equation (2), corresponding to the
determination of the coefficients [a, b, c, d, e, f, h].
The traditional method to determine the coefficients of Equation (2) corresponds to
regression analysis whereas modern techniques of mathematical optimization are developed
over the last decades. The contribution of optimization to geophysics has been interestingly
growing the last decades due to its efficiently in modelling complex natural systems by
determining the best solution from a set of alternative solutions (Goldberg, 1989). But which
is the main advantage of optimization? The answer to this question would be its ability to
reach the optimal solution for any system even under extreme computational environments,
either when data-sets are limited but also when vast data-sets of high diversity require the
determination of a solution describing sufficiently all the samples given. In the following
section the advantages of optimization versus regression would be analysed and the best
solver for optimization process would be selected through a test involving a number of
important deterministic and stochastic algorithms.
3.2 Optimization
In theory, traditional regression analysis -discussed in previous paragraph- is expected to
provide one possible solution for any given equation. Nowadays optimization addresses
the necessity for determining the best solution for data-sets of complex physical systems.
As described previously the effort lies in determining the optimal solution for the
mathematical model of Equation (2), which leads to the implementation of constrained
optimization techniques. The mathematical problem corresponds to the minimization of
misfit represented by the sum of squares of the residuals, between the logarithms of
observed and predicted values (Equation 3). Constrained minimization problems have
some basic pillars which are briefly given as: (1) the existence of a candidate theoretical
solution, to initiate optimization (2) the existence of an objective function, to evaluate
whether minimizing the misfit is achieved (3) a set of linear constraints serving as bounds
for the coefficients’ determination and (4) the determination of convergence criteria
(Rothlauf, 2006).
8 Scientific and Engineering Applications Using MATLAB
It should be noted that for consistency in this example the same objective function, linear
constraints and convergence criteria have been used during the implementation of the
aforementioned solvers.
Techniques used during optimization to exhaust the search space are classified generally in
three classes: (1) Calculus based techniques (2) Guided Random search techniques and (3)
Enumerative techniques (Filho et al., 1994). In this paper calculus based versus guided
random search techniques will be test through comparison of different solvers whereas
enumerative algorithms will be discarded since “they cannot compete to the robustness
race” when compared with the aforementioned techniques mainly due to the characteristics
of their search domains (Said, 2005).
described by researchers (Stoffa & Sen, 1991; Tavakoli & Pezeshk, 2005), making the
aforementioned solvers known for their application in constrained minimization problems
(Goldberg, 1985). GA’s are not limited by restrictive assumptions, concerning the continuity,
the existence of derivatives and the unimodality of the function in terms of computational
geometry. This is an advantage over regression analysis which often determines local
minima as the solution of a given equation (Goldberg, 1989), while at the same time the
suggested solution is highly dependent on a single point of initialization.
In the following paragraph we focus on the initialization, the process of improvement and
the destination of constrained minimization by using stochastic solvers such as GAs. By
keeping the analogy to biological systems a number of chromosomes/strings form the
genetic prescription for the development and the operation of the organism. In our case the
chromosomes/strings are composed by six genes/characters, representing the set of
coefficients of Equation (2).
GAs start with a possible solution or a set of possible solutions corresponding to a theoretic
attenuation curve and continues with optimization in order to determine the optimal
solution for the given data set. In this study both initialization options have been tested,
corresponding either to a single starting point (Simple Genetic Algorithm) or multiple
random-generated starting points (Genetic Algorithm with initial Population Develoment),
forming an initial population of candidate solutions each one satisfying the given linear
constraints for the determination of the coefficients of Equation (2). In order to ensure well-
dispersed and random initial population development, for the adequate representation of
the search space, Latin Hypercube sampling was used (Diaz-Gomez & Hougen, 2006). The
technique was elaborated by Iman et al. (1981) as stratified sampling without replacement,
whereas in recent years risk analysis software employs Latin Hypercube sampling in
preference of Monte Carlo approach for population development.
During optimization every candidate solution/string satisfying the given linear constraints
is evaluated, through the objective function of Equation (3), for its effectiveness in
minimizing the misfit, hence bringing the response value of the system near a desired value.
Within a generation (group of solutions) fitness scaling serves the purpose of ranking each
candidate solution to facilitate the selection of the best solutions that should surviving in the
next generation. In that way, mimicking nature, the fitter solution survives and can be a
parent individual to the next generation whereas worst fit solutions are penalized.
In the present study a number of 200 generations is considered, each one with population
size of 600 individuals/solutions. Stochastic operators like the cross-over, mutation and
survival-of-the-fittest guarantee diversity of the population (Pan, 1995) forcing the GA to
search the solution space intensively, thereby reducing the possibility that the algorithm will
return a local minimum (Goldberg, 1985). In terms of survival, 2 elite individuals/solutions
are guaranteed to survive to the next generation in this study. A crossover fraction equal to
0.8 specifies the percentage of individual/solutions, other than elite children, produced by
crossover in the next generation. Crossover mimics natural recombination between two
parent chromosomes/solutions during which the offspring chromosome/solution has
changed values of specific genes/numbers with respect to the parent genes/numbers. In
this example two point cross-over has been implemented where the selected string is
divided into 3 segments and then 2 segments are exchanged with the corresponding
segments of another string. Mutation, on the other hand, alters a randomly selected
character/coefficient within a string/solution to create a new possible solution. Adaptive
mutation, used in this example, generates new directions in the search space, with respect to
10 Scientific and Engineering Applications Using MATLAB
the last successful or unsuccessful generation, bounded by the linear constraints set for the
coefficients.
By combining the aforementioned stochastic operators, three versions of GAs, aim to test the
relation between diversity of the population and performance, corresponding to
1. Simple Genetic Algorithm (SGA)
2. Genetic Algorithm with initial Population Development (GADP)
3. Hybrid Genetic Algorithm (HGA).
It is noted that Hybrid Genetic Algorithm (HGA) introduces the solution, determined a
priori by a Simple Genetic Algorithm (SGA), to a deterministic solver, which requires the
existence of derivatives, in order to provide local optima once the Genetic Algorithm has
determined the neighborhood of the global optima (Il-Seok et al., 2004). The author included
this enhanced evolutionary algorithm in the comparison to test the efficiency of the
interaction between stochastic and deterministic solvers.
Simulated Annealing is a meta-heuristic algorithm proposed by Kirkpatrick, et al. (1983) and
Cerny (1985) for the determination of global minima. It mimics the physical process where
metals are slowly cooled so that eventually their crystal structure is frozen. The latter state
corresponds to the determination of the optimal solution, using a minimum energy
configuration during its implementation (Bertsimas and Tsitsiklis, 1993).
In this study optimization starts from a randomly generated initial population and a
hypothesis for a parameter, known as Temperature, slowly decreasing from 100° C by a
factor of 0.0059° C in the process of determining the optimal solution. During
implementation each new possible solution is evaluated, for its effectiveness in minimizing
the objective function (Equation 3), then in case of a lower misfit value the suggested
solution is adapted. Simulated Annealing, an important solver in stochastic minimization
problems (Bohachevsky et al., 1986), has been reported to successfully determine global
minima however the author agree with the results of Ingber (1993) that this algorithm
requires fine tuning to specific problems relative to other solvers.
difference however between Simulated Annealing (SA) and deterministic solvers (NLLSQ,
NLLSQDP, PS) is that the latter depend on the existence of derivatives in order to continue
their iterations.
The evaluation of the performance of optimization solvers follows a qualitative and
quantitative method. Analytically by qualitative criteria the authors refers to inability of the
solver to determine a possible solution (1) within the given number of generations (2)
satisfying the linear constraints and (3) in a timely manner. When the above criteria failed,
optimization was implemented again with different starting points, which is acknowledged
to be the main source of error that could lead to a solver’s failure. The alternative solution,
that of relaxing convergence criteria in the case of a solver’s fail, was not considered since it
would jeopardize the final comparison between different solvers. In the event of a solver’s
failure to produce a possible solution for second time, it was excluded from the quantitative
comparison. In that sense Non Linear Least Squares (NLLSQ) solver failing the (2) criterion
has been implemented again using this time multiple starting points (NLLSQDP). After the
second failure to determine a feasible solution by returning the initial conditions -describing
only the theoretical ground motion prediction equation, which was subjectively set by the
programmer- Non Linear Least Squares (NLLSQ, NLLSQDP) have been excluded from the
comparison from this point forward.
Table 1 presents the results of the quantitative comparison of the optimization solvers
together with the coefficients of Equation (2) together with the numerical details used
during optimization, such as the linear constraints introduced in the form of lower and
upper bounds. It is noted that convergence criterion, alternatively known as tolerance, was
set to 1E-06 for the purpose of this test. The quantitative comparison has been based in two
criteria (1) the standard error, in logarithm base 10, and (2) the average sample log-
likelihood (LLH) value (Scherbaum et al., 2009) of the resulting ground motion prediction
equation as derived by a specific solver.
The standard error (σk) has been calculated as the mean of absolute residuals between
observed and predicted by the model gk (where k denotes the index of the solver) ground
motion values described in the equation below:
1
σk = obsi g k ( xi ) (4)
N
Assuming that the set of observations adequately describes nature, the likelihood L(gk|x), of
the model gk given the set of observations x, would represent how close model gk describes
reality. According to Scherbaum et al. (2009) the average sample log-likelihood (LLH)
estimator has been calculated as the mean of log-likelihood values over N number of x
samples
1 N
(
log L ( g k x ) ) = ∑ log ( gk ( xi ) )
N i =1
(5)
It is noted that the sigma value is calculated as the standard deviation of the residuals in log
10 units, for consistency with Equation 2, using the std function. It is noted that the
calculation of log-likelihood estimator of Equation 5 has been made through the function
normlike.
12 Scientific and Engineering Applications Using MATLAB
5. Conclusions
Two major conclusions can be drawn from the results of this study: firstly, deterministic
algorithms, such as Non Linear Least Squares (NLLSQ, NLLSQDP) fail to determine the
whole set of coefficients since the values of the coefficients c, e and h (Table 1) remain fixed
to their lower boundary value. The above remark emphasizes the weakness of deterministic
algorithms leading to the determination of local minima instead of returning a global
solution for the minimization problem. Secondly, the effectiveness of GAs in solving
minimization problems, even in their simpler parameterization (SGA), is supported by their
ranking following the LLH criterion. Thus, the use of Genetic Algorithms aided by initial
Population Development (GADP) by Latin Hypercube sampling for the constrained
minimization problem set in Equation (2) is suggested.
Once the final set of coefficients is determined the seismologist can assess the expected
ground motion at any given site, located at distance R from the epicenter of a strong
earthquake with magnitude M. Figure 8 presents the estimated ground motion after a
magnitude M5 and M6 earthquake at a rock site for various distances using stochastically
derived ground motion prediction equation using Genetic Algorithm with initial Population
Developement.
Ground Motion Estimation During Strong Seismic Events Using Matlab 13
Fig. 8. Attenuation of peak ground acceleration (PGA) over a wide range of distances from
the seismic source.
14 Scientific and Engineering Applications Using MATLAB
Ground motion prediction equations are an important tool for the scientific and engineering
community to forecast anticipated ground motion and predict potential economical losses
and structural damages. Another important application of ground motion prediction
equations lies in developing possible scenarios for the planning short and long term
emergency response.
This chapter describes how Matlab can be used in scientific research for Digital Signal
Processing, Data Archiving but also for modeling complex natural systems through
Optimization. Since the number of graduate students writing computer codes from scratch,
in order to expand the frontiers of their research, continue to grow, core Matlab functions
can be used to develop new software packages in the future. The application examples of
this chapter clearly show that in seismological practice, demanding mathematical
procedures can be implemented and their results can be easily visualized through Matlab.
6. References
Alander, J.T. (2009). An Indexed Bibliography of Genetic Algorithms and Simulated Annealing:
Hybrids and Comparisons, Report for the Department of Electrical Engineering and
Automation, Vaasa, Finland.
Audet, C. & Dennis, J.E. (2003) Analysis of generalized pattern searches. SIAM Journal on
Optimization, Vol. 13, No. 3, pp. 889–903.
Bertsimas, D. & Tsitsiklis, J. (1993). Simulated Annealing, Statistical Review, Vol. 8, No. 1, pp.
10-15.
Bohachevsky, I.O., Johnson, M.E. & Stein, M. L. (1986). Generalized Simulated Annealing for
Function Optimization, Technometrics, Vol. 28, No 3, pp. 209-217.
Cerny, V. (1982). A thermodynamical approach to the travelling salesman problem: An
efficient simulation algorithm, Journal of Optimization Theory and Applications, Vol.
45, pp. 41-51.
Davis, L. ( 1991). Handbook of Genetic Algorithms, Van Nostrand Reinhold, New York.
Goldberg, D. E. (1989). Genetic Algorithms in search, optimization and machine learning,
Addison Wesley Longman, Reading, Massachusetts.
De Jong, K. A. (1975). An analysis of the behavior of a class of genetic adaptive systems.
Dissertation Abstracts International, Vol. 36, No. 10, 5140B.
Diaz-Gomez, P.A. & Hougen, D.F. (2006). Genetic Algorithms for hunting snakes in
hypercubes: Fitness function and open questions, Proceedings of the Seventh ACIS
International Conference on Software Engineering, Aritificial Intelligence, Networking and
Parallel/Distributed Computing, Las Vegas, 2006.
Dolan, E.D., Lewis, R. M. & Torczon, V. (2003). On the local convergence of pattern search,
SIAM Journal on Optimization, Vol. 14, No. 2, 567-583.
El-Mihoub, T.A., Hopgood, A.A.,Nolle L. & Battersby, A. (2006). Hybrid genetic algorithms:
A review, Engineering Letters, Vol. 13, No. 2, pp. 124-137.
Filho, J.L., Treleaven, P.C. & Alippi, C. (1994). Genetic-Algorithm programming
envirinments, Computer, Vol. 27, No. 6, pp. 28-43.
Fogel, L.J., Owens, A.J., & Walsh, M. J. (1966). Artificial intelligence through simulated evolution,
John Wile, New York.
Gabere, M.N. (2007). Simulated Annealing Driven Pattern Search. Algorithms for Global
Optimization, Msc Thesis, School of Computational and Applied Mathematics,
University of the Witwatersrand, Witwatersrand.
Ground Motion Estimation During Strong Seismic Events Using Matlab 15
Grefenstette, J.J. (1986). Optimization of control parameters for genetic algorithms, IEEE
Transactions SMC, Vol. 16, pp. 122-128.
Holland, J. H. (1975). Adaptation in natural and artificial Systems: An introductory analysis with
applications to biology, control, and artificial intelligence, A. Arbor (editor), University
of Michigan Press, MIT press, Cambridge.
Hookes, R. & Jeeves, T.A. (1961). 'Direct search' solution of numerical and statistical
problems, Journal of the Association for Computing Machinery (ACM), Vol. 8, No.2, pp.
212–229.
Il-Seok O., Jin-Seon L. & Byung-Ro M., Hybrid Genetic Algorithms for Feature Selection,
IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 26, No.11, pp.
1424-1437.
Iman, R. L., Helton, J. C., & Campbell, J. E. (1981). An approach to sensitivity analysis of
computer models, Part 1. Introduction, input variable selection and preliminary
variable assessment, Journal of Quality Technology, Vol. 31, No.3,pp. 174-183.
Ingber, L. (1983). Simulated annealing: Practice versus theory, Journal of Mathematical
Computational Modelling, Vol. 18, No.11,pp. 29-57.
Karl, J. H. (1989). An introduction to digital signal processing. Academic Press Inc., San Diego,
California.
Kirkpatrick, S., Gelatt, Jr. C.D, & Vecchi, M.P. (1983). Optimization by simulated annealing,
Science, Vol. 220, No.4598,pp. 671-680.
Kramer, S. L., 1996. Geotechnical Earthquake Engineering. Prentice Hall International Inc.,
New Jersey, pp652.
Liggett, J.A. & Salmon, J.R. (1981). Cubic spline boundary elements, International Journal for
Numerical Methods in Engineering , Vol.17,pp. 543-556.
Mitchell, M. ( 1996). An Introduction to Genetic Algorithms, MIT Press, Cambridge.
Rothlauf, F. (2006). Representations for Genetic and Evolutionary Algorithms, Springer Verlag,
Berlin Heidelberg New York.
Pan, Z., Chen, Y., Kang, L. & Zhang, Y. (1995). Parameter estimation by Genetic Algorithms
For nonlinear regression, Proceedings of International Conference on Optimization
Technique and Applications, G.Z.Liu (editor).
Said, Y. H. (2005). Genetic Algorithms and their applications. Handbook of Statistics 24,
Calyampudi Radhakrishna Rao, Edward J. Wegman, Jeffrey L. Solka (editors),pp.
359-390.
Scherbaum, F., Delavaud, E. & Riggelsen, C. (2009). Model selection in seismic hazard
analysis: An information-theoretic perspective, Bulletin of Seismological Society of
America, Vol.99,pp. 3234–3247.
Segou, M., Voulgaris, N., Makropoulos, K. C., & Stavrakakis, G. N. (2008). A review of the
Greek strong motion database: needs, improvements and future development, in
Proceedings of the 31st general Assembly of the European Seismological Commission,
Hersonissos, Crete, Greece, 7-12 September 2008.
Segou, M. and Voulgaris, N. (2010). Proschema: A Matlab application for processing strong
motion records and estimating earthquake engineering parameters, Computers &
Geosciences, Vol.10,pp. 977-986.
Solomatine, D. P. (1998). Genetic and other global optimization algorithms - comparison and
use in calibration problems, Proceedings of the 3rd International Conference on
Hydroinformatics, Balkema Publishers, Copenhagen, Denmark.
16 Scientific and Engineering Applications Using MATLAB
Stoffa, P. L., & Sen, M. K. (1991). Nonlinear multiparameter optimization using genetic
algorithms: Inversion of plane wave seismograms, Geophysics, Vol. 56, No.11,pp.
1794-1810.
Tavakoli, B., & Pezeshk, S. (2005). Empirical-stochastic ground-motion prediction for eastern
North America, Bulletin of Seismological Society of America, Vol. 95,pp. 2283-2296.
Wetter, M. & Wright, J. (2003). Comparison of generalized pattern search and a genetic
algorithm optimization method, Proceedings of the Eighth International Conference on
Building Simulation (IBPSA), Eindhoven, Netherlands.
2
1. Introduction
Since the earlier times of documented seismological observations, it was noticed that an
earthquake (usually a large one) was followed by a sequence of many smaller earthquakes,
originating in the epicentral region; the first, larger, earthquake is called the mainshock, or
main shock or main event, and the following, smaller, earthquakes are called aftershocks. These
sequences, and their spatial and temporal distributions, depend on the characteristics of the
mainshock and on the physical properties and the state of stress, strain, temperature, etc., of
the region of occurrence (Kisslinger, 1996).
The observation that in many aftershock sequences the magnitude of the largest aftershock
is about ∆M ≈ 1.2 less than that of the mainshock is known as Båth’s law (Helmstetter &
Sornette, 2003; Richter, 1968); it implies an energy ratio E( M aft ) / E( Mmain ) ~ 0.007 between
largest aftershock and mainshock. An earthquake following some mainshock but not small
enough to be considered an aftershock by the Båth’s law criterion is often considered to be a
mainshock in its own right and to constitute, together with the mainshock, a multiple event or
multiplet. A large event followed by an even larger one is demoted from mainshock to
foreshock; a probabilistic calculation, based on the observation that aftershocks follow the
Gutenberg-Richter relation (Gutenberg & Richter, 1954; Kisslinger & Jones, 1991;
Shcherbakov et al., 2005), indicates that for ∆M = 1.2 the seismicity following a given
mainshock has a ~ 6.3% probability of including a “daughter” event larger than the
mainshock (Helmstetter et al., 2003; Holliday et al, 2008).
Thus, the simple definition of aftershock as an earthquake occurring after a mainshock and
in its epicentral region, although implying some causal relation with the mainshock, is
partly semantic and largely circumstantial. Indeed, smallish earthquakes that constitute the
background seismicity occur all the time in a seismic region in the absence of large events, and
continue occurring whether or not a large earthquake occurs, so that not all earthquakes
occurring in the region after a mainshock are necessarily aftershocks.
If aftershocks are a result of the occurrence of the mainshock, then they should be related in
a physical way with its rupture process. The aftershock-producing mechanism is not yet
known, but it is conceivably related with adjustments to the post-mainshock stress field (Lay
& Wallace, 1995) possibly through viscolelastic processes or through fluid flow (Nur &
Booker, 1972); whichever the actual process, aftershocks should be related with the rupture
plane of the mainshock. Kisslinger (1996) qualitatively defines three kinds of aftershocks:
18 Scientific and Engineering Applications Using MATLAB
Class 1 is those occurring on the ruptured area of the fault plane or on a thin band around it.
Class 2 is events that occur on the same fault but outside of the co-seismic ruptured area.
Class 3 is events occurring elsewhere, on faults other than the one ruptured by the
mainshock; these events, whether in the same region or not, will not be considered here as
aftershocks, but rather will be classified as triggered earthquakes.
The number of aftershocks decreases with time after the mainshock according to the
modified Omori relation (Utsu, 1961, 1969, 1970) as
k
N (t ) = , (1)
(t + c ) p
where t is time measured from that of the mainshock, and k, c, and p are positive parameters
which vary with the lithologic, tectonic, and other conditions of the study area. Commonly,
p ranges from 0.9 to 1.8, with most instances between 1.1 and 1.4 (Utsu, 1961); values which
do not show dependence on the magnitude of the mainshock (Utsu, 1962).
Aftershocks occurring within 24 to 48 hours after a large earthquake locate mostly over the
co-seismic rupture area, and provide a good estimation of it (Lay & Wallace, 1995), which
indicates that seismicity at the time is mostly class 1; over longer times the aftershock area
increases (Felzer & Brodsky, 2006; Helmstetter & Sornette, 2002; Mogi, 1968; Tajima y
Kanamori, 1985; Valdés et al., 1982), including at the edges class 2 events.
Among other reasons why aftershock identification is important, we can mention the
following few examples. Aftershocks can give important information about the rupture
area; also, from estimations of co-seismic slip on the fault plane by inversion of seismic
waves, several authors have found that aftershocks are scarce in areas of maximum slip and
concentrate around their edges (Dreger et al., 1994; Engdahl et al, 1989; Hauksson et al., 1994;
Mendoza & Hartzell, 1988), so that aftershocks give information about the rupture process
of the mainshock. Aftershocks can also yield information about the properties of the
epicentral region (Knopoff et al., 1982; Kisslinger, 1996; Kisslinger & Hasegawa, 1991;
Figueroa, 2009), and about possible triggering mechanisms (Roquemore & Simila, 1994).
Since aftershocks can be large enough to contribute to the damage (particularly after
structures have been debilitated by a mainshock), it is important to estimate the hazard
associated with aftershock activity (Felzer et al., 2003; Reasenberg & Jones, 1989).
It has been proposed that, since aftershock activity depends among other factors on the
stress status of the region, there is research on whether some characteristics of the aftershock
activity from intermediate-sized earthquakes can be useful as precursory data for large
earthquake hazard estimation (e.g. Jones, 1994; Keilis-Borok et al., 1980; Wyss, 1986).
Finally, for some studies concerned with large earthquakes, aftershocks can be considered as
noise, and have to be eliminated from the catalogs (e.g. Gardner. & Knopoff, 1974;
Habermann & Wyss, 1984).
In order to use or eliminate aftershocks it is first necessary to identify them. Many methods
have been used, ranging from visual inspection (Molchan & Dmitrieva, 1992) to
sophisticated numerical techniques. A common method identifies aftershocks as those
shocks locating within temporal and spatial windows having lengths which usually depend
on the magnitude of the mainshock (Gardner & Knopoff, 1974; Keilis-Borok et al., 1980;
Knopoff et al., 1982). More sophisticated methods identify aftershocks as belonging to a
spatial cluster, consisting of events within a given distance of at least one other event
belonging to the cluster, which includes the mainshock (Davis & Frohlich 1991a,b; Frohlich
Aftershock Identification Trough Genetic Fault-Plane Fitting 19
& Davis, 1985). A variation of the window method considers events from larger to smaller
magnitudes with the size of the spatial windows a function of the magnitude, and density of
events (Prozorov & Dziewonski, 1982; Prozorov, 1986). Other methods include recognizing
some statistical property (e.g. Kagan & Knopoff 1976, 1981; Vere-Jones & Davies 1966) or
interpreting the relations between events according to some statistical chain or branching
model (Molchan & Dmitrieva, 1992; Reasenberg, 1985).
Our method includes some of the above mentioned techniques used to discard events which
cannot be aftershocks, and then proceeds to identify aftershocks based on the physical
model of a rupture plane and on recognized statistical relationships. An early
unsophisticated application of the rupture plane model, which proved that this principle of
aftershock identification was feasible, was part of an unpublished MSc. thesis (Granados,
2000).
2. The method
We work with seismic catalogs containing occurrence time (days), hypocentral x (East), y
(North), and z (up), and, optionally, horizontal and vertical location uncertainties uh and
uv , all these in kilometers. If location uncertainties are not included in the catalog, optional
horizontal and vertical uncertainties are assigned equally to all events.
Any events occurring before the event with the largest magnitude Mmax , the mainshock, are
eliminated. All spatial coordinates are then referred to those of the mainshock.
A rough time cutoff, eliminates events occurring after more than an optional cutoff time
(default is nyr = 4 years ), because after a few years it is difficult to distinguish aftershock
activity from background seismicity.
The extent of the aftershock area depends on the energy, i.e. on the magnitude, of the
mainshock (Utsu, & Seki, 1955), as does the rupture area. A first rough spatial
discrimination, based on an average of the empirical magnitude M vs. source-length r (km)
relationships of Wells and Coppersmith (1994):
or (Kagan, 2004):
r = 20 ⋅ 10 ( M − 6)/2 ; (2b)
eliminates events farther away from the hypocenter than 1.5 times the r f length estimated
by (2).
Next, a spatial clustering analysis, where events separated by no more than a given critical
distance r, of the order of hundreds of meters to a few kilometers, depending on the spatial
coverage of the catalog, are considered to be related, is used to eliminate events which do
not relate to the mainshock or to other possible aftershocks.
The parameters in the modified Omori’s law (1) are not known a priori; they are estimated
from the statistics of the aftershocks (Davis & Frohlich, 1991b; Guo & Ogata, 1997; Ogata,
1983), which we do not yet have. However, this relation tells us that for long enough times
after the occurrence of the mainshock the number of aftershocks decreases until seismic
activity in the epicentral region returns to its background level (Ogata & Shimazaki, 1984).
20 Scientific and Engineering Applications Using MATLAB
When aftershock occurrence shows gaps comparable to those characterizing the background
seismicity, we can consider that the aftershock activity is, if not ended, at least scarce
enough to be comparable to the background activity and can no longer be distinguished
from it. The critical gap length depends on the region and the magnitude threshold of the
observations; we use a default critical gap length t gap = 10days . The gap is measured as the
average of a given number of inter-event times (default ngap = 10 ). Events occurring after a
critical gap are discarded from the possible aftershocks.
Weights, based on relation (1), are optionally assigned to the remaining events, using typical
values for c (default c ≈ 2 ) and p (default p ≈ 1.0 ) and by setting k = c p ⇒ N (t = 0) = 1 , so
that those events which have a large likelihood of being aftershocks have weights ~1, while
later events have smaller weights:
cp
wi = . (3)
(ti + c ) p
Next, plane fitting is carried out iteratively; at each iteration, a plane that passes through the
mainshock hypocenter is fitted to all remaining events, through a genetic scheme described
below, and fit outliers (events too far away from the plane) are discarded. Iteration
continues until the goodness-of-fit criterion is met (successful fit) or until a preset maximum
number of iterations is attained (unsuccessful fit).
The t gap , c, and p values can be refined using the final results of a first, tentative aftershock
determination, to do a second one.
The genetic plane search for the plane, characterized by its azimuth φ and dip δ , which
minimizes the L1 norm of the perpendicular distances from the fault plane to all events, is as
follows. An initial set of equispaced φ j and δ j values, covering the acceptable ranges, is
built and the fit to the aftershock candidates is evaluated for each pair of values.
To estimate the error of fit for each candidate plane, for each azimuth φ j and dip δ j pair,
event coordinates ( x , y , z ) are rotated as:
and in the ( x′, y′, z′) coordinate system the equation of the plane is x′ = 0 , so that the L1 fit
error for the j’th parameter pair is easily computed as
Na Na
ε j ≡ ε( φ j , δ j ) = ∑|x′i (φ j , δ j ) wi | / ∑ wi , (5)
i =1 i =1
Next, N p ( N p − 1) / 2 children are constructed with parameters which are averages of those
of each pair of parents; other N c children are constructed for each parent j by randomly
modifying one or the other of the parameter values, according to
(
φ j N p + k = N φ j , σm
φ ε j / εNp ); k = 1,… , N c , (6)
δ j Np + k = N( δ ,σj
m
δ εj / εNp )
{ }
φ = max σmin , σφ , σδ = max {σ min , σδ } , and
where N designates the normal distribution, σm m
where f σ is a damping factor (default is 1.25) and ui is the location uncertainty, adjusted
for the plane dip as
N ( t ) = ∑ H ( t − t i ) K i ( t − t i + c i )− p ,
i
22 Scientific and Engineering Applications Using MATLAB
t
N ( t ) = ∑ χ i K i ( t − t i + c )− p , (9)
i t k − ti
where
t 0≤t≤1
χ(t ) = Π(t − 1 / 2) = ,
0 other t
Π(t ) is the boxcar function, and tk > ti is the time of the succeeding event occurring along
the same fault plane and cluster or the time when the gap criterion is met ( tk = ∞ if no such
event or gap occur). An event at t k may “inherit” some of the aftershocks from the one at
ti , but they will be identified as aftershocks and be duly eliminated. Båth’s law indicates
that ∆M should be 1.2, but this criterion results, for most catalogs, in too many mainshocks;
our default value is thus ∆M = 1.0 , which means that we accept as aftershocks those with
energy ratio E( M ) / E( Mmain ) < ~ 0.016 , but this ∆M value can be easily changed by the
user. The largest events in the catalog are plotted vs. time above the initial and resulting
event densities, to illustrate which aftershocks were eliminated. The program iterates the
whole process, as many times as needed, until no more aftershocks are found
In cleancat, parameters are not set interactively, but can be easily adjusted in a list of
adjustable parameters at the beginning of the code.
The program optionally outputs a catalog excluding identified aftershocks.
3. Application
We will now show some examples of the application of the method. The aftplane program
will be used to identify fault planes and aftershocks from three mainshock-aftershock
sequences from two different parts of the world featuring different tectonic environments.
The cleancat program will be used to clean the catalog of a fault system.
parameters c = 2.0 days and p = 1.0 ; aftershock magnitude criterion ∆M = 1.0 ; fit criterion
maximum error εmax km.
Fig. 1. Seismicity map of southern California showing the location of the Joshua Tree and
Landers faults (both within the red diamond.) (Modified from Lin et al, 2007.)
Fig. 2. Joshua Tree mainshock plus 3497 acceptable aftershock candidates (left) and Joshua
Tree mainshock plus 3379 clustered aftershock candidates (right).
24 Scientific and Engineering Applications Using MATLAB
Fig. 3. Joshua Tree mainshock (red asterisk) plus 1094 aftershocks (blue diamonds); left:
plan view showing 171.6° faultplane azimuth; right: view along faultplane azimuth clearly
showing the 86.6° dip.
The main Joshua Tree shock and the identified 1094 aftershocks are shown in figure 3, both
in a plan view (left) which clearly shows the resulting 171.6° faultplane azimuth, and a cross
section along the fault plane azimuth (right) which shows the resulting 86.6° faultplane dip.
The values found by aftplane agree extremely well with those estimated by Velasco et al.
(1994) of strike 171°, dip 89°. Figure 4 shows a cross section parallel to the fault plane,
illustrating aftershock concentrations.
Fig. 4. Joshua Tree mainshock (red asterisk) plus 1094 aftershocks (diamonds), cross section
seen along azimuth 81.6°, perpendicular to the 171.6° faultplane azimuth.
Aftershock Identification Trough Genetic Fault-Plane Fitting 25
Fig. 5. Landers mainshock plus 17553 acceptable aftershock candidates (left) and Landers
mainshock plus 12834 clustered aftershock candidates (right).
The main Landers shock and the identified 3,225 aftershocks are shown in figure 6, in a cross
section seen along the determined 340.6° fault plane azimuth, which shows the resulting
70.1° faultplane dip. The values found by aftplane agree extremely well with those
estimated by Velasco et al. (1994) of strike 341°°, dip 70°. Figure 7 shows a cross section
parallel to the fault plane, illustrating aftershock concentrations.
Fig. 6. Landers mainshock (red asterisk) plus 3225 aftershocks (blue diamonds), view along
340.6° faultplane azimuth. The location of the mainshock hypocenter is obscured by those of
the aftershocks. Dip 70.1°
Fig. 7. Landers mainshock (red asterisk) plus 3225 aftershocks (blue diamonds), view along
azimuth 70°, perpendicular to 340.6° faultplane azimuth.
26 Scientific and Engineering Applications Using MATLAB
Fig. 8. Location of the Armería, 22 January 2003, MW 7.3 earthquake; the star indicates the
epicenter of the main shock, circles are located events following the mainshock, located by
the RESCO network.
The main Armería shock and the identified 460 aftershocks are shown in figure 10, in a
cross section seen along the determined 86.2° fault plane azimuth, which shows the
resulting 33.2° faultplane dip. The values found by aftplane agree extremely well with the
above mentioned estimates strike 300° and 20° to 30° (Nuñez et al., 2004; Yagi et al., 2004;
Andrews et al., 2010). Figure 11 shows a cross section parallel to the fault plane, illustrating
aftershock concentrations.
Aftershock Identification Trough Genetic Fault-Plane Fitting 27
Fig. 9. Armería mainshock plus 10868 acceptable aftershock candidates shown as black
crosses (left) and Armería mainshock plus 7109 clustered aftershock candidates shown as
blue circles (right).
Fig. 10. Armería mainshock (red asterisk) plus 460 aftershocks (blue diamonds), view along
86.2° faultplane azimuth.
Fig. 11. Armería mainshock (red asterisk) plus 460 aftershocks (blue diamonds), view along
azimuth 356.2°, perpendicular to 86.2° faultplane azimuth.
28 Scientific and Engineering Applications Using MATLAB
Fig. 11. Joshua Tree- Landers fault system seismicity (black crosses) and identified
aftershocks (yellow circles).
Figure 12 shows the occurrence times and magnitudes of the largest events in the catalog
(top), and below them (middle) is plotted the ocurrence density (per ∆t = 46 day) vs. time
(blue line); peaks in the occurrence rate after large shocks aftershocks are clearly seen. The
bottom panel of figure 12 shows the occurrence density after aftershock elimination using
∆M = 1.0 (red line); although densities are much lower than before filtering, the peaks
coinciding with the occurrence of the Joshua Tree and Landers earthquakes indicate that
many aftershocks are not being identified.
Use of ∆M = 0.9 (Fig.13) effectively diminishes the troublesome above mentioned peaks to
background seismicity level; the seismicity rate peak around t ~ 7.25 10 5 days is associated
with distributed seismicity with events about the same size. Thus, we see that Båth’s law is
not appropriate for the mainshock-aftershock relationships in the seismicity of the Joshua
Tree-Landers fault system; a small (0.1 unit) change in the magnitude criterion can make a
large difference in the aftershock recognition capability.
Aftershock Identification Trough Genetic Fault-Plane Fitting 29
Fig. 12. Joshua Tree-Landers fault system seismicity versus time. The top panel shows the
largest events in the period (blue circles with vertical lines). The middle and bottom panels
show seismicity rates, for 46 day-long time intervals, before (middle) and after (bottom)
processing by cleancat, respectively; note the different vertical scales.
Fig. 13. Joshua Tree-Landers fault system seismicity rates, for 46 day-long time intervals,
after processing by cleancat, with ∆M = 0.9 .
4. Conclusions
We present a simple method for identification and/or elimination of aftershocks, based on
the generally accepted assumption that aftershocks are related to the fault rupture of the
mainshock. The method has been tried on various catalogs with good results and, when
aftershocks are numerous enough, good estimates of rupture planes that agree very well
with those reported in the literature.
A variation of the method used for eliminating all aftershocks from a seismicity catalog
(“catalog cleaning”) uses, iteratively, a variation of the same principle. Using the seismicity
occurrence time rate as illustration and criterion of the effectiveness of the method, indicates
that the required difference between mainshock and aftershocks ∆M is a key parameter for
correct aftershock identification, and that ∆M = 1.2 (Båth’s law) may be too strict for some
geographic areas and/or seismo-tectonic settings.
5. Acknowledgements
Many thanks to José Frez, Juan García A., and María Luisa Argote for useful criticism and
comments. We are grateful to RESCO and Gabriel Reyes, and to the SCEC for the use of
their catalogs. We also thank Mr Davor Vidic of Intech for the kind invitation to participate
in the present book.
30 Scientific and Engineering Applications Using MATLAB
6. References
Andrews,V., Stock, J. Ramírez-Vázquez,C., & Reyes-Dávila, G. (2010). Double-difference
Relocation of the Aftershocks of the Tecomán, Colima, Mexico Earthquake of 22
January 2003. Pageoph, DOI 10.1007/s00024-010-0203-0.
Davis, S. D. & Frohlich, C. (1991a). Single-link cluster analysis, synthetic earthquakes
catalogues and aftershock identification. Geophys. J. Int. 104, 289-306.
Davis, S. D. & Frohlich, C. (1991b). Single-link cluster analysis of earthquake aftershock:
decay laws and regional variations. J. Geophys. Res. 96, 6335-6350.
Dieterich, J. (1994). A constitutive law for rate of earthquake production and its application
to earthquake clustering. J. Geophys. Res. 99, 2601-2618.
Dreger, D., Pasyanos, M., Loper, S., McKenzie, R., Gregor, N., Uhrhammer, B., &
Romanowicz, B. (1994). Source process of the 17 January 1994 Northdridge
earthquake. EOS, 75(16), 103. (abstract)
Engdahl, E. R., Billington, S., & Kisslinger, C. (1989). Telesismically recorded seismicity
before and after the May 7, 1986, Andreanof Islands, Alaska, earthquake”. J.
Geophys. Res. 94, 15481-15498.
Felzer,K., Abercrombie,R., & Ekstrom,G. (2003). Secondary aftershocks and their importance
for aftershock forecasting. Bull. Seismol. Soc. Am. 93, 1433-1448.
Felzer,K., Abercrombie,R., & Ekstrom,G. (2004) A common origin for aftershocks,
foreshocks, and multiplets. Bull. Seismol. Soc. Am. 94, 88-98.
Felzer,K. & Brodsky,E. (2006). Decay of aftershock density with distance indicates triggering
by dynamic stress. Nature 441, 735-738.
Figueroa, A. (2009). Analysis of inter-event time in aftershock sequences for identification
stress relaxation status., MSc Thesis, UNAM, 89pp (in Spanish).
Frohlich, C. & Davis, S. D. (1985). Identification of aftershocks of deep earthquakes by a new
ratios method. Geophys. Res. Lett. 12, 713-716.
Gardner, J. & Knopoff, L. (1974). Is the sequence of earthquakes in Southern California with
aftershocks removed Poissonian? Yes. Bull. Seismol. Soc. Am. 64, 1363- 1367.
Granados,J. (2000). Identification of seismic aftershocks in catalogs. MSc.Thesis in Earth Sciences,
CICESE, 145pp (in Spanish).
Guo, Z. & Ogata, Y. (1997). Statistical relations between the parameters of aftershocks in
time, space, and magnitude. J.Geiphys.Res.102(B2), 2857-2873.
Gutenberg, B. & Richter, C. (1954). Seismicity of the Earth and associated phenomena. Princeton
Univ. Press, 310pp.
Habermann, R. and Wyss, M. (1984). Background seismicity rates and precursory seismic
quiescence: Imperial Valley, California, Bull. Seismol. Soc. Am. 74, 1743–1755.
Hauksson, E., Hutton, K., & Kanamori, H. (1994). The Mw 6.7 Northdridge, California,
earthquake of January 17, 1994 and its aftershocks. Program Northdridge Abstr., 89th
Annu. Meet. Seismol. Soc. Am., unpublished abstract.
Helmstetter, A. & Sornette, D. (2002). Diffusion of epicenters of earthquake aftershocks,
Omori’s law and generalized continuous-time random walk models. Phys.Rev.E 66,
061104.
Helmstetter, A. & Sornette, D. (2003). Båth’s law derived from the Gutenberg-Richter law
and from aftershock properties. Geophys.Res.Lett.30, 2069.
Aftershock Identification Trough Genetic Fault-Plane Fitting 31
Helmstetter, A., Sornette, D., & Grasso, J. (2003). Mainshocks are aftershocks of conditional
foreshocks: How do foreshock statistical properties emerge from aftershock laws.
J.Geophys.Res. 108(B1), 2046.
Holliday, J., Turcotte, D., & Rundle, J. (2008). A review of earthquake statistics: fault and
seismicity-based models, ETAS and BASS. Pageoph 165, 1003-1024.
Kagan, Y. & Knopoff, L. (1976). Statistical search for non-random features of the seismicity
of strong earthquakes. Phys. Earth planet. Inter. 12, 291-318.
Kagan,Y. (2004) Short-Term Properties of Earthquake Catalogs and Models of Earthquake
Source. Bull. Seismol. Soc. Am. 94, 1207–1228
Keilis-Borok, V., Knopoff, L., & Rowain, I. (1980). Bursts of aftershocks long term precursors
of strong earthquakes. Nature. 283, 259- 263 p.
Kisslinger, C. (1996). Aftershocks and fault-zone properties. Advances in geophysics, 38, 1-36.
Kisslinger, C. & Hasegawa, A. (1991). Seismotectonics of intermediate-depth earthquakes
from properties of aftershocks sequences. Tectonophysics 197, 27-40.
Kisslinger, C. & Jones, L. M. (1991). Properties of aftershocks sequences in southern
California. J. Geophys. Res., 96, 11947-11958.
Knopoff, L., Kagan, Y., & Knopoff, R. (1982). b-values for foreshocks and aftershocks in real
and simulated earthquake sequences. Bull. Seismol. Soc. Am. 72, 1663-1675.
Jones, L. M. (1994). Foreshocks, aftershocks, and earthquake probabilities: Accounting for
the Landers earthquake. Bull. Seismol. Soc. Am., 84, 892 – 899 p.
Lay, T. & Wallace, T. (1995). Modern global seismology. Academic Press, Inc., 521pp.
Lin, G., Shearer, P. M., & Hauksson, E. (2007). Applying a three-dimensional velocity
model, waveform cross correlation, and cluster analysis to locate southern
California seismicity from 1981 to 2005. J. Geophys. Res. 112, B12309,
doi:10.1029/2007JB004986.
Lomnitz, C. (1966). Magnitude stability in earthquake sequences. Bull. Seismol. Soc. Am. 56,
247- 249.
Mendoza, C. & Hartzell, S. H., (1988). Aftershock patterns and mainshock faulting. Bull.
Seismol. Soc. Am. 78, 1438-1449.
Mogi, K. (1968). Development of aftershock areas of great earthquakes. Bull. Earthquake Res.
Insti., Tokyo Univ., 46, 175 – 203.
Molchan, G. M. & Dmitrieva, O. E. (1992). Aftershock identification: Methods and new
approaches. Geophys. J. Int. 109, 501-516.
Núñez-Cornú, F. J., Reyes-Dávila, G. A., Rutz-López, M., Trejo-Gómez, E., Camarena-
García, M. A., & Ramírez-Vazquez, C. A. (2004). The 2003 Armería, México
Earthquake (M w 7.4): Mainshock and Early Aftershocks, Seism. Res. Lett. 75, 734–
743.
Nur, A. & Booker, J. (1972). Aftershocks caused by pore fluid flow? Science, 175, 885-887.
Nur,A., Ron,H. & Beroza,G. (1993). The nature of the Landers-Mojave earthquake line.
Science 261, 201-203.
Ogata, Y. & Shimazaki, K. (1984). Transition from aftershock to normal activity: the 1965 Rat
Islands earthquake aftershock sequence. Bull. Seismol. Soc. Am. 74, 1757-1765.
Ogata, Y. 1983. Estimation of the parameters in the modified Omori formula for aftershock
frequencies by the maximum likelihood procedure. J. Phys. Earth, 31, 115 –124.
32 Scientific and Engineering Applications Using MATLAB
1. Introduction
Every year during the summer months (June-September), the southern part of Asia, in
particular the Indian subcontinent receives continuous and widespread rains due to the
monsoon (meaning seasonal change of wind direction) more specifically known as the
summer or South-West (SW) monsoon. The phenomena in summer takes place due to the
cross hemispheric reversal of winds bringing in considerable amount of water vapour from
the high pressure regions over the relatively colder Indian and Pacific oceans and the
Arabian sea to the low pressure system over heated land mass areas. Due to large
geographical coverage and high inter annual variability in terms of associated rainfall it
constitutes one of the important elements of the global climate system (Mooley &
Parthasarathy, 1984). Also the overall mean monsoon precipitation distribution significantly
depends on the intra-seasonal oscillations of dry and wet spells (Waliser et al., 2003). The
sudden onset of monsoon is characterised by a highly energised pattern of lightning,
thunderstorm, cloud burst and incessant rainfall over a large area of southwest coastal
region of the Indian Kerala state (Soman & Kumar, 1993). Normally the first episodic rain of
the monsoon occurs over Burma and Thailand in the middle of May and then extends to the
northwest over most of India within a month. The northward progression of monsoon is
symptomatic of a large-scale transition of deep convection system (called the tropical
convection zone) from the oceanic-equatorial to tropical-continental regions (Sikka & Gadgil
1980).
While the SW monsoon period (June-September), accounting for about 80% of the total
rainfall is vital to meet the Indian agricultural and hydrological requirements, the winter
monsoon flowing from North East direction during October-November contributes only
marginally albeit meeting the ground water requirement of some areas falling in the shadow
region of the summer monsoon cover. Many oceanic and atmospheric parameters like the El
Nino and Southern Oscillation (ENSO), Eurasian snow cover in winter/spring, northern
hemispheric temperature in winter etc. influence the year-to-year variability of the monsoon
rainfall over India (Rajeevan et al., 2004). Hence the inter annual variations severely affect
the ground water resources of not only India and neighbouring countries like Sri Lanka,
Bangladesh and Pakistan but also on small scales, the equatorial Africa, northern Australia,
and south-western United States.
34 Scientific and Engineering Applications Using MATLAB
The monsoon rainfall intensity has intra- and inter- annual variation at different spatial
scales impacting on the availability of water for agriculture and other uses. As a result
almost half of the world's population living in monsoon region suffer from food and water
insecurity. The detailed scientific reasons of such variation in monsoon rainfall conditions
are not clearly understood. Effort directed to make accurate prediction of the overall
monsoon rainfall of the season averaged over the whole country as such is quite useful for
farming activity. More desirable objective of such forecasts is to predict the prospects of
seasonal rainfall with finer details of its spatial distribution so that the peasants may decide
about the sowing activity at district or even village levels for better crop yields.
There are three major monsoon variables: (i) the yearly onset date of monsoon over India (it
first enters the coastal Kerala state by 1st June ± 8 days and later progresses to cover the
whole country), (ii) total seasonal rainfall during June-September (the area weighted
summer monsoon long term average rainfall for the whole country as estimated by India
Meteorological Department (IMD) is about 88 cm with a coefficient of variation of 10%
which in real terms may lead to widespread drought/flood) and (iii) seasonal rainfall
distribution over different Indian meteorological sub-divisions (prediction of this is the most
difficult). To gain an early knowledge of the amount of seasonal monsoon rainfall or
determination of its trend has been a challenging scientific problem for long. A number of
empirical, statistical and dynamical/general circulation models have been developed for the
purpose by various groups the world over (Munot & Kumar, 2007). For the Indian region
mainly the statistical models have been used on an operational basis with partial success.
For example the statistical model could not forecast the recent rain deficient years of 2002,
2004 and 2009. Moreover the variations of the 2 monsoon indices namely the time of onset
and the total amount of monsoon rains are not directly correlated to the rainfall distribution
in different parts of the country which would follow the dynamics of meteorological
parameters at local or micro levels.
The Asian–Australian monsoon (another name of the same Indian summer monsoon in a
regional context) is a coupled geophysical phenomenon the intensity of which is regulated
through negative feedbacks between the land, ocean, and atmosphere. Indian Ocean
heat transport calculations using 41-yr (1958–98) data revealed that the Indian Ocean heat
transport possesses strong variability at all time scales from intra seasonal (10–90 days) to
inter annual (a biennial signal is significant). The amplitude of the intra seasonal variability
is similar to the seasonal cycle, and the amplitude of the inter- annual variability is about
one-tenth of the seasonal cycle (Chirokova & Webster, 2006).
According to IMD’s definition of monsoon transition, at the surface the onset is recognized
as a rapid, substantial, and sustained increase in rainfall over a large spatial scale while the
withdrawal marks the return to dry, quiescent conditions. The criterion is that the rainfall
amounts over Kerala district increase from below 5 to over 15 mm/day during the onset
(Anathakrishnan and Soman, 1988). A different condition that assesses the onset and
withdrawal dates of the Indian monsoon has also been derived from variability in the large-
scale hydrologic cycle. The hydrologic cycle is chosen as a key physical basis for monitoring
the monsoon due to the essential roles played by zonal and meridional gradients in water
vapour, clouds, and rainfall in driving the large-scale monsoon circulation. Lateral
transports of water vapour are required for the sustenance of monsoon rains. To diagnose
onset and withdrawal, vertically integrated moisture transport (VIMT) is considered more
representative by some authors instead of rainfall, which over the large scale is often poorly
Sea Surface Temperature (SST) and the Indian Summer Monsoon 35
measured and modelled. An index, named the hydrologic onset and withdrawal index
(HOWI), is thus formed from those regions where VIMT variability is pronounced at
the beginning and end of the monsoon season respectively. Analysis of inter annual
variability in monsoon onset and withdrawal dates based on the HOWI reveals
robust associations that are weak and insignificant when assessed using other onset criteria.
The HOWI criterion shows strong correlations between total rainfall and both onset and
withdrawal of monsoon (Fasullo & Webster, 2006). But these indices have a drawback of
being determined retrospectively and hence their less predictive potentials.
While the tropical day-to-day weather conditions have a restricted predictability of 2–3
days, the seasonal (June-September) mean monsoon circulation in the tropics is
potentially more predictable (Rajeevan, 2001). This is understandable as the low frequency
or longer period oscillatory features of the tropical variability is caused by slowly varying
boundary conditions and forcings due to sea surface temperature (SST), land surface
temperature, soil moisture, snow cover, etc. (Charney & Shukla, 1981). But a considerable
fraction of monsoon variability results from internal dynamics at higher frequencies (intra
seasonal) often due to local/regional effects of environmental factors including aerosol
distribution, atmospheric pollution, orography, forest dynamics etc. This intra seasonal
variation is the main limiting factor of the monsoon predictability at subdivision level. Still,
overall it makes sense in estimating and forecasting the likely onset date of the summer
monsoon as well as the total seasonal rainfall in three categories of normal (~area weighted
average value of ~88 cm, excess (~97 cm) and deficient (~79 cm). This information is very
important for macro level planning of ‘Kharif’ crop cultivation with paddy as the main crop.
The ‘Kharif’ crops such as paddy, sugarcane, groundnut, maize, pulses etc., need timely and
adequate water either through rains or through artificial irrigation system. Indian
agriculture still depends heavily on monsoon rains and sowing times differ with locations
and with crop-type during April-July months. It would therefore be ideal to get early alerts
at micro level, a difficult proposition at present but may be realized in future. In absence of
this, the accurate predictions of the monsoon onset dates (date over Kerala governs onsets
over other regions following a climatological pattern of monsoon progress as shown in Fig-1
for the south Asian region and over India) and the integrated seasonal rainfall, help in
managing the agricultural output to a large extent.
Efforts for accurately predicting the onset date and the overall strength of seasonal
precipitation have been continuing for a long time using synoptic data analysis, empirical,
statistical and dynamical modelling by the IMD (Hastenrath & Greischar, 1993; Raghu Kant
&. Iyengar, 2003). Out of a number of parameters SST anomalies of the Pacific/Indian ocean
related to the strong El Nino and La-Nina events and associated circulation like ENSO have
been correlated with delayed/weak and early/strong monsoon rainfall over India (Joseph et
al, 1994; Nakazawa, 1988; Philander, 1990).
Due to the availability of all-weather and homogenized SST data sets since 2002-03 from
microwave sensor (AMSRE) of AQUA satellite, it has been possible to carry out a
quantitative assessment of SST anomalies with 3-day temporal and 0.25ºx0.25º lat-long grid
(pixel unit) resolutions. Hence a more detailed investigation of the effect of SST on the
monsoon can be tested both for variations in the onset dates and also the total seasonal
rainfall over monsoon fed regions. As a demonstration to utilizing the satellite data for
monitoring the SST vis-à-vis the monsoon system over India, a study has been carried out to
develop a real time model and the preliminary results published (Chakravarty, 2009).
36 Scientific and Engineering Applications Using MATLAB
Fig. 1. Climatological dates of the onset of the south Asian summer monsoon (adapted from
Fasullo and Webster, 2002). The monsoon onset date contours are also shown over India.
In addition to SST data from satellites, upper air balloon and rocket experiments have been
regularly conducted from Thiruvananthapuram, the capital of the Kerala State to detect any
changes in the temperature and wind fields up to ~70 km owing to monsoon circulation. A
special campaign called ROMEX (Rocket Monsoon Experiment), involving balloon and
rocket borne wind measurements was carried out during April-June, 2007 for real-time
monitoring of middle atmospheric zonal/meridional winds to study the prognostic
potential of early reversal of upper winds for the possible date of setting in of summer
monsoon over the Kerala coastal region (Chakravarty & Namboodiri, 2007). The essential
statistics for such prediction was built up using voluminous data already collected between
1971-90 using balloon/rocket flights from Thumba Equatorial Rocket Launching Station
(TERLS) of the Indian Space Research Organisation (ISRO), Thiruvanathapuram. The results
showed distinct potential of upper tropospheric and stratospheric wind reversal parameters
(circulation indices) being used for predicting the monsoon onset day from the
climatological trends of these circulation indices and 2007 campaign data used as a test case.
In the background of the above summary providing the present status of the various
observational and modelling aspects and their applications to predict the prospects of
seasonal monsoon rainfall, the main purpose of this chapter is to (a) review and carry out
further studies on the development of a real time model using SST as the main parameter
over a region covering central/western Pacific ocean and the Indian ocean to predict at
regular intervals the prospects of possible onset date and total rainfall over India (b) as a
collateral and supporting information to (a) above, demonstrate possible use of upper air
parameters like the troposphere/stratosphere circulation indices for predicting the monsoon
onset and strength for the ensuing season, (c) possible operationalisation of the MATLAB
programme with features to provide periodic updates on the monsoon rainfall during April-
September with built in graphics for various plots to show the monsoon trends. This
assumes that the AQUA/AMSRE satellite type SST and other data would continue to be
available.
SOI often indicate El Niño episodes. These negative values are usually accompanied by
sustained warming of the central and eastern tropical Pacific Ocean, a decrease in the
strength of the Pacific Trade Winds, and a reduction in rainfall over eastern and northern
Australia. Though there are departures as the strong El Niño event of 1997/98 had only
marginal effect on Australia. But severe droughts resulted from the weak to moderate El
Niño events of 2002/03 and 2006/07. Positive values of the SOI are associated with stronger
Pacific trade winds and warmer sea temperatures to the north of Australia, popularly
known as a La Niña episode. Waters in the central and eastern tropical Pacific Ocean
become cooler during this time. Together these give an increased probability that eastern
and northern Australia will be wetter than normal. The strong and preceding El Niño/La
Niña episodes are also used mainly to parameterise a statistical correlation coefficient. This
aspect is further dealt in this chapter. The method used by the Australian Bureau of
Meteorology is based on the Mean Sea Level Pressure (MSLP) difference between Tahiti and
Darwin. It is calculated as follows:
[ Pdiff − Pdiffav ]
SOI = 10
SD( Pdiff )
where Pdiff = (average Tahiti MSLP for the month) - (average Darwin MSLP for the month),
Pdiffav=long term average of Pdiff for the month in question, and SD(Pdiff) = long term
standard deviation of Pdiff for the month in question. The multiplication by 10 is a
convention. Using this convention, the SOI ranges from about –35 to about +35, and the
value of the SOI can be quoted as a whole number. The SOI is usually computed on a
monthly basis, with values over longer periods such as year being sometimes used. Daily or
weekly values of the SOI do not convey much in the way of useful information about the
current state of the climate and can fluctuate markedly because of daily weather patterns,
and should not be used for climate purposes. The monthly average SOI data are used to
calculate the annual and pre-monsoon period (Jan-April) mean SOI values for 1970-2010
period. The monthly average SOI contours are generated with respect to the different years
to identify El Nino and La-Nina episodes. The official monsoon onset dates over Kerala for
different years are taken from daily weather bulletins, Met. Center, Thiruvananthapuram,
Kerala for correlation studies with ENSO, SST etc.
The global daily SST data from AQUA satellite’s TMI/AMSRE sensors are produced by
Remote Sensing Systems and sponsored by the NASA Earth Science REASoN DISCOVER
Project and the AMSR-E Science Team. Data since 2002-03 are available at their official
website (www.remss.com). The daily global SST values for the period 2003-11 are examined
during the pre-monsoon as well as monsoon months. A MATLAB computer program is
developed to read daily binary AMSRE files of global SST (3 days aggregate values at a
fixed morning time every day) and generate global contour plots. After going through a
large number of day’s data of all these years a part of Pacific Ocean region bounded by -20°
to 10° in latitude & 80° to 240° in longitude is selected as test site for further analysis. This
oceanic region (ignoring a small part of land area) is termed Nino-Broad Pacific or Nino-BP
to distinguish it from other existing Nino region definitions. Fig-2 shows the areas bounded
by existing Nino regions and the new Nino-BP region used in the present analysis. The
program then can focus on any Nino area and produce time series of SST pixel (25 km x
25km area) number values distributed in the selected 27-31 ºC temperature bins with a
38 Scientific and Engineering Applications Using MATLAB
resolution of 1 ºC separately for the years 2003-11. Total number of pixels within this test
region is 65,000. Using MATLAB graphics, the SST contours and bar charts (at weekly
intervals) are automatically plotted separately for each year covering Jan-May and June-
September months. The programme with the regular input of TMI/AMSRE data can thus be
used to follow the trends for prospects of monsoon strength on real time basis.
Fig. 2. Oceanic regions called Nino-1, Nino-2, Nino-3, Nino-4, Nino-3, 4 for carrying out
different sensitivity analysis and the new Nino-BP region defined here for association with
Indian monsoon sensitivity.
3. Results
The results obtained are categorised in the following sub sections: (a) variations of annual
average, pre-monsoon and monthly averages of SOI and its correlation with onset and
strength of monsoon rainfall, (b) global and Nino-BP region maps of SST, (c) variations of
pixel-based structure of SST in the Nino-BP and Nino-4 regions and efficacy of a real time
model for monitoring the monsoon system and (d) upper air indices of monsoon circulation.
SOI smoothened by taking 12 months running means of the monthly average SOI data for
the same period. A high value of correlation coefficient (~0.84) is found between the annual
average and pre monsoon months average values indicating that the effect of these events
during pre monsoon months normally continue in a similar manner of warming or cooling
for a longer period of the year. Fig-4 shows the contour plot of the monthly mean ENSO
index values for the period. It is seen that the El Nino and La Nina events take place in an
alternating manner both having periodicities of ~5-6 years. The darkest (black) shades
indicate occurrence of El Nino and lightest (white) shades La Nina episodes. Other grey
shades indicate different intermediate states of these events including the normal condition
of neither being present. There are higher frequency seasonal structures within the episodic
year that can be seen in the contours.
Fig. 3. Variation of ENSO index during the period 1970-2010. Top panel shows the annual
and pre monsoon months (Mar-May) average values and the bottom panel 12 months
running means.
The determination of the monsoon onset date over Kerala is normally announced and later
slightly modified based on the event definition by IMD. Since there is another criterion of
the hydrologic onset and withdrawal index (HOWI) mentioned earlier, a comparison is
made by plotting these two sets of dates in Fig-5. The HOWI index data is available only
during 1970-2000. The figure shows some major differences between the onset dates with a
correlation coefficient of ~0.61. Hence the studies related to the causes of variations of onset
dates would have less meaningful results if there were a change in the definition itself. It is
found that a burst of rainfall due to non-monsoon reasons may happen locally to mimic an
early onset of monsoon applicable to both the data sets with differing degrees.
40 Scientific and Engineering Applications Using MATLAB
Fig. 4. Contour plot of monthly mean ENSO index during 1970-2010. The darker shades
show El Nino and the lighter La Nina episodes.
Notwithstanding the imperfect way the monsoon onset is defined, it may still be helpful to
examine any direct influence of the El Nino/La Nina events on the monsoon onset dates
over Kerala coast. The annual and pre-monsoon period means of SOI are plotted with
respect to the onset dates of each year during 1970-2010 in Fig-6. From the distribution it is
Fig. 5. Comparison of monsoon onset days using IMD and HOWI data sets.
Sea Surface Temperature (SST) and the Indian Summer Monsoon 41
Fig. 6. Scatter plot of mean values of SOI and monsoon onset dates for the years 1970-2010.
clear that only in 4 cases out of 40 years of data, the SOI exceeded 15 in the negative scale
(strong El-Nino event) and for all of these years the monsoon onset dates were delayed, i.e.,
beyond the normal onset date of 1 June. For all other values of SOI ranging between –14 and
+15 the correlation is very poor. There have been much more delayed onset for smaller
negative values of mean SOI. It may be noted that at least for positive values of SOI the
onset has taken place within 4 June. So from normal conditions to La-Nina events provide a
better prognosis that the monsoon would be timely compared to the El Nino events
association with delayed monsoon. Apart from the onset date the other parameter is the %
deviations of total area weighted seasonal (June-September) rainfall from long term mean
value over India. Table-1 shows a few selected years of strong El Nino or La Nina events
and associated SOI indices, onset dates and % deviations of rainfall. It is noted that there is a
better correlation of negative and positive association to rainfall index with strong El Nino
and La Nina events respectively.
From the above analysis it is noted that the mere presence of El Nino/La Nina events does
not in itself provide any definite clue to the prospects of the ensuing monsoon (onset date in
particular). However the normal conditions of oceanic state and mild La Nina situations
appear to be favourable for a better monsoon. It is not possible to quantify these effects on
the basis of these events. Hence it is necessary to deal with the oceanic state in a more
quantitative manner. The sea surface temperature variations can be studied in more details
using regular satellite data at a high resolution. The main parameter which is responsible for
rainfall under monsoon system is the evaporation from the ocean surface which is closely
linked to cloud formation and rainfall while being transported away from the source region.
Fig. 7. Global SST contour map using AQUA data On 9 April 2007. Blue to red colour
coding is to distinguish low and high temperature pixels. There are no measurements over
the land areas and these are shown in white. Concentration of high temperature pixels in
the Indian and western Pacific oceans is a nominal phenomena for April month.
Sea Surface Temperature (SST) and the Indian Summer Monsoon 43
The same programme can be employed to focus on the specific Nino regions of interest and
map the SST values for different days. As an illustration for the Nano-BP region the SST
maps are drawn on 3 days in the pre monsoon months (Mar-May) for two years of 2004 and
2008. Fig-8 shows these 2 sets of SST contours for comparison. The year 2004 had an early
monsoon onset on 18 May and year 2008 a near normal onset on 29 May. The total seasonal
rainfall over India for 2004 was 12% deficient and 2008 almost normal. It can be noted from
the figure that the SSTs were generally higher during 2004 compared to 2008 for the same
days of the year. The figure also shows that there is a larger region in the eastern Pacific
Ocean which is colder during 2008 with a gradual decrease in the size of this area as the
days progress compared to those during 2004. More quantification of these results is
required and this has been attempted in terms of the distribution of SST values at pixel
levels.
Fig. 8. SST contours of 3 selected days during pre-monsoon period for the years 2004 and
2008 over the Nino-BP region. Red/blue colours indicate cold/hot regions. White areas are
landmasses.
well behaved from the view point of overall monsoon seasonal rainfall distribution over
India. From Fig-9 it is clear that the basic trend of the curves are similar but there is a large
variation in absolute values of mean SST from year to year. The lowest temperatures were
observed in 2008 (a near perfect monsoon year) and the values up to April for the current
year (2011) also shows low SSTs indicating that the monsoon during 2011 could be like 2008
unless there are some drastic systemic deviations during the monsoon months.
Fig. 9. Mean SST values over Nino-BP region for different days of pre-monsoon months at
weekly interval for individual years (2004-11).
Fig. 10. Mean SST values over Nino-BP region for different days of monsoon months at
weekly interval for individual years (2004-10).
Sea Surface Temperature (SST) and the Indian Summer Monsoon 45
For other 4 years which include 2 normal years (2007, 2010) and 2 abnormal years (2004,
2009) no definite conclusions can be drawn based on the trend and absolute values of mean
SST. Fig-10 provides another clue that mean SSTs should come down to 2008 level for the
monsoon to pick up. It is noted that while the mean SSTs showed gradual cooling up to end
August for 2007 and 2010 the cooling trend was not enough for 2004 (12% deficient rainfall
occurred compared to long period average, LPA) and 2009. In fact for 2009 mean SSTs
remained high resulting in a major failure of monsoon rains (23 % deficient from LPA).
Hence the mean SST over Niano-BP region can provide a better handle by monitoring it on a
real time basis which may help in short-term (a week or so in advance) prediction.
In order to check if the actual pixel level temperatures can provide additional information
the contours of number of pixels distributed over the temperature bins 25-31 °C for the
observation days (at weekly interval) for the pre monsoon period for 4 selected years (2004,
2007, 2008, 2009) are drawn and shown in Fig-11. The pattern of variation of number of low
and high temperature pixels are similar in pairs for [2007, 2008] and for [2004, 2009]. Both
2004 and 2009 years show presence of high number of pixels ≥ 28 °C. This would have
implications in losing the water vapour flux transport through oceanic and pre-monsoon
rains before the onset of monsoon. Sometimes this may mimic arrival of monsoon at an
early date. This appears to have happened during 2004.
Fig. 11. Contour plots of number of pixels having different temperature values between 25-
31 °C during pre-monsoon months for different years, 2004, 2007, 2008, and 2009 over Nino-
BP region.
46 Scientific and Engineering Applications Using MATLAB
Fig. 12. Mean SST values over Nino-4 region for different days of pre-monsoon months at
weekly interval for individual years (2004-11).
2004 2007
2008 2009
Fig. 13. Contour plots of number of pixels having different temperature values between 25-
31 °C during pre-monsoon months for different years, 2004, 2007, 2008, and 2009 over Nino-
4 region.
Similar analysis has been carried out for the Nino-4 region and the results presented in Fig-
12 and Fig-13. The trends of mean SST over Nino-4 are similar but being a small region it
could get influenced by local effects. From Fig-13 the similarities in paired years of [2007,
2008] and [2004, 2009] cannot be easily discerned as in the case of Nino-BP plots. However it
Sea Surface Temperature (SST) and the Indian Summer Monsoon 47
would help to conduct SST analysis for both Nino-BP and Nino-4 regions as part of the
monsoon forecasting and real time modelling.
For determining the characteristics of SST pixel distribution during the build up phase of the
monsoon system, monthly average values are computed and shown as bar charts over Nino-
BP for the years 2004, 2007, 2008 and 2009 in Fig-14. While excess number of pixels with
SST≥ 28 °C is found for anomalous years 2004 and 2009 particularly in May, these were
found to be only little lower during 2007. Similar analysis in Fig-15 for the Nino-4 area
shows a clearer demarcation of pixel distribution comparison between the two pairs of
years. During 2010 there was a strong El Nino effect from Jan-Mar and hence the prospects
of monsoon looked bleak. But the situation improved subsequently with a La-Nina setting
in. This change over took place at right time so that monsoon rains picked up from July
onwards. A comparison of relevant parameters between 2008 and 2010 is shown in Fig-16. It
can be noted that the pixel number distribution with SST in both the years are similar for the
month of May which changed the prospects for a better monsoon during 2010.
2004 2007
2008 2009
Fig. 14. Average monthly pixel (25-31 °C) number distribution for pre monsoon months of
Jan-May for 2004, 2007, 2008 and 2009 over the Nino-BP region.
2004 2007
2008 2009
Fig. 15. Average monthly pixel (25-31 °C) distribution for pre monsoon months of Jan-May
for 2004, 2007, 2008 and 2009 over the Nino-4 region.
48 Scientific and Engineering Applications Using MATLAB
2008 2010
Fig. 16. Comparison of pixel number contour and monthly average pixel distribution over
Nino-BP region during 2008 and 2010.
Fig. 17. (Clockwise) Comparison of monthly mean Nino-BP region pixel distribution for
2008 and up to April of 2011, time series of mean SST over Nino-BP region for 2008 and up
to April, 2011 and contours of pixel numbers pertaining to temperature bins of 25-31 °C up
to April of 2011.
Sea Surface Temperature (SST) and the Indian Summer Monsoon 49
Fig-17 shows a comparison of monsoon build up parameters of 2008 and 2011 (data used up
to April 2011). It clearly shows a striking similarity in the variations of the mean monthly
pixel distribution, the absolute values and the trend of mean SST over Nino-BP region
between 2008 and 2011 so far (April, 2011). The contour of pixel numbers for SST bins also
shows similarity to the 2008 plot up to end of April. Hence all conditions are favourable for
2011 monsoon to be a success and there are chances that it would be like the one during
2008 along with its near uniform spatial distribution.
Fig. 18. Vector wind profiles of individual days of observation using RH-200 rocket and
rawinsonde balloon launches from Thumba during April-June 2007 under ROMEX
campaign.
Fig. 19. Smoothed contour plot of zonal winds between 0-70 km height ranges over Thumba.
Red contours show westerly winds and blue easterlies.
for such strengthening. Acquiring a value of mean zonal wind of 10 m/s in this
troposphere/ lower stratosphere height range appears to be a rough figure to work with for
linking with monsoon.
Sea Surface Temperature (SST) and the Indian Summer Monsoon 51
Fig. 20. Time series of mean zonal winds for height ranges, 11-20 km, 41-50 km and 51-60 km
during ROMEX campaign period.
Fig. 21. Trends of mean zonal winds for 11-20 km and 41-50 km height ranges from long
period (1971-90) balloon/rocket data differentiated for normal (blue) and delayed (red)
monsoon and how the ROMEX data for 2007 fares as a test case (black).
Fig. 21 shows the statistical data of mean zonal winds between 11-20 km and 41-50 km
distinguished in terms of normal and delayed monsoon with red and blue lines respectively.
Over these climatological pattern of variation the ROMEX line for 2007 is superimposed. It
is clear that the 2007 data fits closer to the blue lines which are for the normal monsoon
onset days.
52 Scientific and Engineering Applications Using MATLAB
4. Conclusion
4.1 The rainfall during the SW monsoon period (June-September) constitutes the main
source of water in India for agriculture (particularly pertaining to the 'Kharif' crops),
hydroelectric power generation and drinking water requirements. There is a large inter
annual and intra seasonal variation of this water resource due to variations in the monsoon
genesis and progress. While the spatial distribution of monsoon rainfall is caused due to fast
response parameters and local impacts, both the onset date and the season's integrated
rainfall are caused by slowly varying boundary conditions and forcings like that of the sea
surface temperature, snow cover etc. Hence the variation of sea surface temperature is taken
up here for prognosis studies.
4.2 Based on the AQUA satellites all weather microwave sensor data a near real-time
interactive computer model has been developed to extract daily minimum global SST values
of 1440x720 pixels, each pixel covering 25 km x 25 km of lat-long area. The programme also
selects specific oceanic region like the Nino-4 and Nino-BP to compute daily mean SST over
the region and can add into previous days data to generate a real time trend. Such trends of
previous years during 2003-2010 are used to study the variations and its influence on the
onset dates and the seasonal rainfall. In the background of the statistics or the real-time
model the data of current year is analysed up to April 2011. The progress of the absolute SST
values and its trend indicate a near normal monsoon during 2011 somewhat similar in
characteristics that of 2008.
4.3 There have been many studies to link the occurrence of major El Nino and La Nina
events with the changes in the monsoon onset date. Such analysis is carried out in this
report also and it is found that the strong El Nino/La Nina events have a negative/positive
effect on monsoon. But considering the whole range of seasonal mean SOI it is found that
statistically the positive SOI values have positive impact but negative SOI values have both
positive and negative impacts. Thus the SOI index has only limited applicability as a
predictor parameter. Hence the main purpose of this study is to make a more quantitative
assessment at pixel level of the SST linked monsoon variability.
4.4 A novel sensitivity analysis is carried out by selecting a broader Pacific Ocean region
called Nino-BP (defined by the author). Within this region the pixels are counted and placed
in temperature bins of values 25-31 ºC. The resultant matrix provides daily number of pixels
distributed over these temperature bins. This pixels numbers are plotted as time-bin
contours or as bar graphs. The main result from this analysis shows that larger number of
pixels distributed in temperature bins ≤ 28 ºC during pre monsoon or during monsoon
months has a positive impact on the onset date and total seasonal rainfall. Reverse applies
for larger number of pixels in temperature bins of value > 28 ºC.
4.5 As collateral and useful information on the prospects of onset date, the time sequence of
mean zonal wind values over TERLS between 11-20 km and 41-50 km height ranges provide
the possible transition date about 4-5 days in advance by comparing its real time trend with
a statistical model of long period balloon and rocket data during 1971-90.
5. Acknowledgement
AQUA/AMSR-E data are produced by Remote Sensing Systems and sponsored by the
NASA Earth Science REASoN DISCOVER Project and the AMSR-E Science Team. Data
are available at http://www.remss.com. Author acknowledges the help provided by Dr. K.
Sea Surface Temperature (SST) and the Indian Summer Monsoon 53
6. References
Ananthakrishnan, R & Soman, M. K. (1988). The onset of southwest monsoon over Kerala,
1901–1980. J. Climatol., Vol.8, pp. 283–296
Chakravarty, S. C. (2009). A computer model to study the variability of grid-based Sea
Surface Temperature (SST) values derived from AQUA/AMSRE satellite data and
its influence on the onset of South West Monsoon near the Kerala coastal region in
India. IEEE Xplore, pp. 1-5, ISBN: 978-1-4244-4562-2
Chakravarty, S. C.; Datta, J. & Revankar, C. P. (1992). Climatology of long-period oscillations
in the equatorial middle atmosphere over Thumba, India. Current Science, Vol.63,
No.1, pp. 33-42
Chakravarty, S. C. & Namboodiri, K. V. S. (2007). A Scientific Report on the Rocket
Monsoon Experiment (ROMEX) campaign (2007). ISRO, India (unpublished)
Charney, J. G. & Shukla, J. (1981). In: Monsoon Dynamics, J. Lighthill, (Ed.), Cambridge
University Press, Cambridge, pp. 99–110.
Chirokova, G. & P. J. Webster (2006). Interannual Variability of Indian Ocean Heat
Transport. Journal of Climate, 19, 1013-1031
Fasullo, J. & Webster, P. J. (2002). A hydrological definition of Indian Monsoon Onset and
withdrawl. Journal of Climate, Vol.16, pp. 3200-3211
Fukao, S. (2006). Coupling processes in the equatorial atmosphere (CPEA): A project
overview. J. Meteor. Soc. Japan, Vol.84A, pp. 1-18
Hastenrath, S. & Greischar, L. (1993). Changing predictability of Indian monsoon rainfall
anomalies? Proc. Indian Acad. Sci. (Earth Planet. Sci.), Vol.102, pp. 35–47
Joseph, P. V.; Eisheid, J. & Pyle, R. J. (1994). Interannual variability of the onset of the Indian
summer monsoon and its association with atmospheric features, El Nin˜o, and sea
surface temperature anomalies. J. Climate, Vol.7, pp. 81–105, 1994
Mooley, D. A. & Parthasarathy, B. (1984). Fluctuations in all India summer monsoon rainfall
during 1871–1978. Climatic Change, Vol.6, pp. 287–301
Munot, A. A. & Krishna Kumar, K. (2007). Long Range prediction of Indian summer
monsoon rainfall. J. Earth Sys. Sci., Vol.116, No.1, pp. 73-79
Nakazawa, T. (1988). Tropical super clusters within intraseasonal variations over the
western Pacific. J. Meteor. Soc. Japan, Vol.66, pp. 823–839
Philander, S. G. H. (1990) In: El Niño, La Niña and the Southern Oscillation, Academic Press,
San Diego, CA, pp. 1-289
Raghu Kanth, S. T. G. & Iyengar, R. N. (2003). Empirical modeling and forecasting of
Indian monsoon rainfall. Current Science, Vol.85, pp. 1189-1201
Rajeevan, M. (2001). Prediction of Indian summer monsoon: Status, problems and prospects.
Current Science, Vol.81, No.11, pp. 1451-1457
Rajeevan, M. D.; Pai, S., Dikshit S. K. & Kelkar, R. R. (2004). IMD’s new operational models
for long-range forecast of southwest monsoon rainfall over India and their
verification for 2003. Current Science, Vol.86, No.3, pp. 422-431
54 Scientific and Engineering Applications Using MATLAB
Sikka, D. R. & Gadgil, S. (1980). On the maximum cloud zone and the ITCZ over Indian
longitudes during the southwest monsoon. Mon. Wea. Rev., Vol.108, pp. 1840–1853
Soman, M. K. & Kumar, K. K. (1993). Space–time evolution of meteorological features
associated with the onset of Indian summer monsoon. Mon. Wea. Rev., Vol.121, pp.
1177–1194
Waliser, D. E.; Stern, W. F., Schubert, S. D. & Lau, K. M. (2003). Dynamic predictability of
intraseasonal variability associated with the Asian summer monsoon. Q. J. R.
Meteorol. Soc., Vol.129, pp. 2897-2925
4
1. Introduction
The Ob River occupies a central place in Western Siberia. Its basin comprises all spatial and
dynamic variability typical for Western Siberian ecosystems. The river floodplains are the
most dynamic parts and in the same time host most of the human activities. However, at
present spatial planning of economic activities in the floodplain areas cannot be based on
water regime and geomorphologic processes, because relevant inventory of relief and
hydrological monitoring data are hardly accessible. Moreover, for nature conservation, in
terms of planning, area selection and management, and data of biodiversity of floodplain
ecosystems are missing.
This research intends to fill these gaps in knowledge. In addition factual relations between
ecosystem productivity of different floodplain units and hydrological regime need to be
studied. There to a hydro-ecological monitoring scheme for analyses of the Ob-Irtysh
floodplain of the region Middle Ob will be set up.
The research will be based on the already developed hydro-ecological zone maps and
hydrological data of hydro-meteorological stations (river stages dynamics) representative
for different river floodplain sections.
The monitoring results of a period of more 50 years stage observations flood hydrographs
have been plotted and frequencies and duration of floods and relations with
geomorphologic characteristics of floodplain relief and flood depths been calculated. These
data may need to be evaluated more thoroughly and thereafter compared with ecosystem
productivity data, which will be gathered within the scope of this project. The bird
population is considered as the most sensitive element of the ecosystem.
Birds are an important component of ecosystems. They function as the consumers of the first
and second orders in the trophic chain of an ecosystem. The main factor determining the
annual dynamics of the bird population in the West Siberia is migration. In autumn the
majority of bird population leaves for south, in spring they come back, and the time of
return coincides with the period of spring high water in the Ob. This period also coincides
with the breeding stage. The floodplain of the Ob attracts birds in the first place. It is related
to the warming action of the waters that the river brings from south to north and to the
more productive and diverse habitat conditions than in the interfluve territory. The Ob
valley serves as a kind of air channel, along which the majority of birds moves. The years
56 Scientific and Engineering Applications Using MATLAB
when the water-content parameters are close to the mean annual values are especially
favorable, while those with low or extremely high water-content parameters are unfavorable
(Adam A. M. & Bolotnov, 1982-2010). Thus, from the point of view of the monitoring of the
state of the floodplain ecosystem the birds serve as a good indicator of its state in terms of
hydro-thermal conditions. In addition to that, they occupy an important place in the nature
management of the region, since a significant part of the waterbird population is a major
hunting resource actively used by the local people and hunters from other regions.
Spring level increase normally starts in the second part of April (early and late periods are
the beginning or the end of April), even while freezing over. In general water level flow of
the Ob river we observe one wave with intensive increase and very slow reduce. In the area
of Tom river inflow we can observe crested water flow or two-three weaked waves of flood,
appeared as a result of split of multy-peaked flood in the upper reaches of the Ob river.
Duration of the floodplain can last from 120 days while “friendly” springs till 150 days. An
average duration of a low meadow flood can consists of 8-12 days by Kruglikovo
hydrometric station, 63 days by Kolpashevo hydrometric station and 68 days by
Alexandrovskoye hydrometric station. The longest period of middle Ob meadow flood lasts
for 2-3,5 months. Although in the north from Alexandrovskoye water point given water
horizons can be observed for 2-2,5 months later than around Kruglicovo hydrometric station
in a period of long springs. The end of floodplain normally comes in July or August. An
average height of water level increase above pre-flood period is equal to 5vm (before Tom
river inflow), than 7-8 and than the highest 9-11 (Alexandrovskoye hydrometric station).
Duration of level increase period is about 30-35 days with average insensitivity of increase
as of 30-35cm per 24 hours.
Analysis of yearly water point level shows that the whole meadow of the Ob river is filled
with water when floodplain reaches the maximal level 10-1% of provision. That’s why we
rarely see the common flood, it happens once per 30-50 years or even more seldom. At the
same time the low meadow of the Ob river from Kolpashevo hydrometric station to the
northern boundary goes under water every year. So the low part of the meadow is more
adopted for floodplain impact and low water floods or floods for short periods can be born
badly. This situation is quiet rare as the Ob river (in the lower area of the Ket river inflow)
exists in a natural regime and high side inflow, huge water cumulative basin, which forms
water reserves while autumn-winter season. These reserves provide obligatory spring
floodplain. When analyzing the flood the meadow area is usually divided into high, middle
and low areas. This dividing is equal to levels which are higher than 25 % of provision and
average multy-year frequency of flooding once per 4 years, close to 50 % of provision and
flooding once per 2 years, and lower than 75 % of provision, with every year flooding. Often
this division happens not objectively – by vegetation which is a vivid indicator of flooding.
However the usage of counting characteristics of provision lets apply fixed marks of the
flood level and refer them to meadow relief and find not only the fact of flooding and also
give quantity characteristics of its duration, height of filling and starting and finishing date.
It is very important for the low meadow areas which are flooded every year but also have
differences in parameters. Dividing of the low meadow territory within the boundaries of
100–95 % of provision permits basically perform the part from the whole annually flooded
territory.
For this purpose we used the long-term observations (1977-2000) conducted in spring and
summer, when the influence of spring high waters and other ecological factors on the
spatial-temporal structure of the bird population was studied in detail. The average bird
population in the floodplain varies from 1000 (willow forests) to 47 ind/km2 (river). The
highest population density is recorded in the villages on the river banks, which varies from
1500 in the first half of summer to 4000 ind/km2 in the second half. The value of the
parameters decreases as the complexity of habitats diminishes, the relief lowers, and the
moisture level grows (from forests and shrubs to the meadows of high ridges, meadows of
depressions, lakes and watercourses). In the second half of summer almost in all habitats the
bird abundance increases 1.5 times. In the forest habitats the yellow-breasted bunting, coal
58 Scientific and Engineering Applications Using MATLAB
tit, and long-tailed tit dominate in abundance, and in the shrub habitats, the reed bunting.
The yellow-breasted and reed buntings dominate in meadows, too. In the over-wet and wet
meadows and lakes the Pallas's grasshopper warbler and garganey are abundant. The sand
martin dominates on the Ob and outlets. In the villages the Eurasian tree sparrow and barn
swallow dominate. In total about 128 bird species live in the area, which can be divided into
6 ecological groups by their habitats: forest-shrub birds, 59 species; birds of dry meadows, 9
species (the common quail, skylark); birds of wet meadows, 12 species (the corn crake, great
snipe, common snipe); water-bog birds on the over-wet and flooded meadows, 8 species
(the Eurasian bittern, mallard); birds of the waterbodies, 21 species (the common teal,
common pochard, tufted duck); birds associated with villages, 11 species (the barn swallow,
Eurasian tree and house sparrows) (Vartapetov, 1984, Yudkin, 1987). The group of birds of
wet meadows is the most dynamic by the value of the year-to-year changes.
stock influences another through the flows. Two main loops that influence the bird
population density value are represented in Figure 1. The upper loop determines the flow of
density increase (FDI), the lower, the flow of density decrease (FDD). NFDI and NFDD are
normal flow of density increase and normal flow of density decrease corresponding to the
mean annual conditions.
Fig. 3. Structure of the model of the dynamics of the bird population density in the
floodplain of the middle Ob.
interrelations between the variables (factors) included in the model. The following
parameters are chosen as the stocks forming the system structure: bird population density
(P) and anthropogenic factors (AF). The marks of irregular shape ("cloud-like")—inflows or
outflows—are positioned outside the system. Any closed loop is a feedback loop. The stock
introduced to the system (AF) reflects the rational human activity that should lead to a
positive effect, therefore, AF provokes the growth of bird population. The share of the
human interference with nature that results in negative consequences (poaching, nest
devastation, change in the natural habitats) leads to a decrease in the bird population. These
phenomena are marked with the "predation" variable (PR). There is a feedback between the
stocks (P) and (AF), too. Its idea is that with high bird population density the number of
The Analysis of Influence of River Floods on
Biotic Components of Floodplain Ecosystems with the Help of MATLAB Simulation 61
birds affected by human activity grows. The feedback is realized with the multiplier (MFPI),
which increases the flow of the AF increase (FAFI) or leaves it unchanged depending on the
density population value. The change of the stocks of the spring high waters by human
activity is expressed with the multiplier of dependence of the stock on the anthropogenic
factor (MSAF).
through the multipliers MIMG and MDMG respectively. The area of the flood-free territory
determines such vital conditions as the presence of territory for nesting, trophic resources
and death by predating. With high flood the birds concentrate on the flood-free area (FA),
and the waterbirds inflow. The share of the chicks dying of predation increases, and the
reproduction success is lowered by the overpopulating. Thus, there are dependencies of FDI
and FDD on FA. In the model they are realized through the multiplier of dependence of FDI
on FA (MFIF) and multiplier of dependence of FDD on FA (MFDF) (Figure 3 shows only
some graphs of dependencies).
The FDI and FDD are influenced by the high water duration (HD), which is a function of the
water level (WL). Its value modifies the FDI and FDD through the multipliers of
dependence of FDI and FDD on HD (MIHD, MDHD) (see Fig. 3c,d).
A sharp drop in temperature (frosts) in the nesting period leads to the decrease in bird
population due to the death of chicks and clutches. To consider the impact of temperature
we introduced the variable of temperature regime (TR) into the model, which influenced the
FDD through the multiplier of dependence of FDD on temperature (MFDT) (Fig. 3e).
One of the factors limiting the bird population is predation, which is especially dramatic
during the nesting period. Under predation we mean an immediate effect on the birds of the
preying animals and an indirect human influence that promotes it (depriving the nests their
defense devices during hay-making, disclosing the nests and hatches by troubling, etc.). In
the model predation is represented with the (PR) variable, which influences the FDD
through the multiplier of dependence of FDD on PR (MDPR) (Fig. 3f). Predation depends to
a certain degree on the flood-free area (FA), and, as was mentioned above, increases as the
areas suitable for nesting decrease. This relation is expressed with the multiplier of
dependence of PR on FA (MPRFA).
The model structure includes an auxiliary variable, phytomass of meadows (PM), which
expresses the height of meadow vegetation and occupied area (Shepeleva, 1986). If its values
are low or high, the bird population decreases. The parameter of phytomass of meadows
(PM) influences the FDD through the multiplier of dependence of FDD on phytomass
(MDPM), and on FDI, through the multiplier of dependence of FDI on phytomass (MIPM).
It is known that the meadow productivity in the floodplain is influences by the duration of
the flooding of the floodplain during high water [10]. In the model it is expressed through
the multiplier of dependence of the phytomass increase on the high water duration (MPHD)
(Fig. 3k).
Migrations conditioned by the character of spring high waters are typical of the birds
inhabiting the floodplain. The larger is the number of birds claiming a nesting area, the
greater is the influence of the spring high waters.
Duration and high level of the high waters result in an increase in waterbird population
accompanied by the general tendency of decrease in the bird population density. Very low
high waters decrease a share of water and near-water birds and cause migration of the dry-
meadow species from the interfluve part into the floodplain. The optimal state of the bird
population is observed in the years with high but short high waters. In the model the bird
migration is represented as the variable (MG), and dependencies are expressed through the
multipliers of FDI and FDD on migration (MIMG and MDMG).
approximation. The outside variables of the model PR, PM, TR, MG, WL are defined as the
functions of time t. The population density of birds in any moment of time is defined as the
density in the antecedent moment of time plus the density added due to FDI and minus the
density decreasing due to FDD in the embraced period.
where F FDI t,t + 1 , is the flow of the density increase on the following interval, ind./km2;
Pt , population density in the given moment, ind./km2; NFDI, normal flow of the density
increase, 1/t; MIAF, the multiplier of dependence of the flow of increase on anthropogenic
factors; MIFA, the multiplier of dependence of the flow of increase on the flood-free area;
MIHD, the multiplier of dependence of the flow of increase on the high-water duration;
MIMG, the multiplier of dependence of the flow of increase on migrations; MIPM, the
multiplier of dependence of the flow increase on the phytomass.
The flow of the density decrease is a part of the reversed feedback. The basic flow of
decrease equals the population density P, multiplied by the normal flow of the density
decrease NFDD. The real flow of the decrease depends on the conditions in the other parts
of the system. Amelioration, predation, hydrological and temperature regimes, and
phytomass of meadows influence the FDD with the multipliers. The equation of FDD is as
follows:
FDDt,t + 1 Pt NFDD • MFDA • MDHD • MFDT • MDPR • MDPM • MFDA • MDMG, (3)
where FDDt,t + 1 is the flow of the density decrease on the following interval, ind./km2;
MFDA, the multiplier of dependence of the flow of decrease on the flood-free area; MDHD,
the multiplier of dependence of the flow of decrease on the high-water duration; MFDT, the
multiplier of dependence of the flow of decrease on temperature; MDPR, the multiplier of
dependence of the flow of decrease on predation; MDPM, the multiplier of dependence of
the flow of decrease on phytomass; MFD A, the multiplier of dependence of the flow of
64 Scientific and Engineering Applications Using MATLAB
where AF t , AF t -1 are the value of the anthropogenic factors in the present and previous
moments of time; FAFI t -1, t , the flow of the anthropogenic factor increase on the previous
interval 1/t.
The flow of anthropogenic factor increase FAFI equals the basic flow or, in this case, normal
flow of anthropogenic factor increase multiplied by the multiplier of dependence of FAFI on
population density (MFPI). This multiplier in normal conditions is equal to 1 and begins to
work in extreme situation, when the population density of birds drops sharply:
where FAFI t , t -1 is the flow of the anthropogenic factor increase on the following interval
1/t; NFAFI, the normal value of the anthropogenic factor 1/t; MFPI t , the value of the
multiplier MFPI in the present moment of time.
Let us consider the mathematical description of the auxiliary variables: variable WL is the
function of time and is given a priori: WL = F(t), variable TR is also a predictable function of
time TR = Ф(t).
In the model the variable PR is determined by its value in normal conditions and state of
two multipliers in the given moment of time, i.e.
where PRt and NPR are the parameters of predation in the present moment of time and
corresponding to the normal
conditions; MPRFAt , MPRAF t , the multipliers of dependence of predation on the flood-
free area and anthropogenic factor.
The duration of the high-water and the flood-free area depend only on the water level:
HD ( WL ), FA ( WL ), MG ( WL ) , (7)
The variable A in the model is represented by the relative value Sa/S, where Sa is the area of
the ameliorated lands, S, the area of Kolpashevskii raion.
The parameter of phytomass in the present moment of time is determined by its normal
value multiplied by the multiplier of dependence of the meadow phytomass on the high-
water duration (MPHD):
MPHDt , the multiplier of dependence of the phytomass on the high-water duration in the
given moment of time.
The equation of the initial conditions is written in the following form: t0 is the initial
reference point; P0, initial density of bird population, ind./km2.
The values of FDI, FDD and FAFI necessary for the first calculation of the model are as
follows:
Fig. 6. Basic model of the dynamics of the bird population in the floodplain of the Ob
realized with the help of MATLAB 5.2.1.
Pto Po ; (12)
AF to AFo ; (16)
FDDt,t + 1 Pt NFDD • MFDA • MDHD • MFDT • MDPR • MDPM • MFDA • MFMG ; (19)
WL F(t ); (22)
HD ( WL ); (23)
FA ( WL ); (24)
TR (t ); (25)
MG ( WL ); (26)
A f(t ); (28)
The quantitative presentation of the stock, flow and auxiliary variables is based on the
experimental data on the real system. When determining the constants and variables, the
conditions of 1977 were taken as the reference points, i.e., the state of the system is described
as related to this year. The dynamics of the population density of birds is followed for 1977-
2000, with spring-summer period considered within each year conditions. The step of
modeling is accepted as equal to one year. All variables of the model are characterized with
relative values.
Fig. 10. Ratio of basic values characterizing the state of an ornithocomplex, as anthropogenic
factor changes: 1, AF; 2, WL; 3, P.
3. Conclusion
The partial change in hydrological regime and industrial use of the flood-plain lands with
moderate amelioration does not affect the dynamics of birds. If the area of amelioration
grows to 50% of the total area of the flood-plain, the existing ecosystem will be destroyed in
the lowest high-waters (75% of provision and less) and will not be able to restore for 4 years.
It is necessary to control the scale of amelioration and not to allow the system to begin
irrevocable destruction.
The built model is of theoretical and applied character. The structure of the model can be
used as basic for biotic components of the flood-plain ecosystem when predicting the basic
tendencies of their behavior and monitoring. It is built for the component which plays an
indication role. Introduction of certain changes into the parameters of water regime can help
in determining the upper and lower limits. When they are passed, the flood-plain ecosystem
begins to change in general.
70 Scientific and Engineering Applications Using MATLAB
For conservation of ecosystem values in the Ob river floodplain the following aspects
should considered:
- Preferably land use types with low impact should be developed: (eco)-tourism,
recreation, trade, small scale agriculture (diary, pastures)
- Preservation an equally balanced land use between natural and semi-natural
ecosystems, given the ecological potentials
- Hay-making is most suitable land use for wet meadows
- Land reclamation development should focus on high floodplain parts
- Development of health-improving recreation.
The conservation of natural resources is achieved by combination of two units. The first unit
provides annual observation of high water regime of Middle Ob flood-land in comparison
with long-term data. Second unit is human activity management, which includes the
preparation of recommendations for the main resource users: administrators, farmers,
hunters and fisherman. It is expected that by the management measures to be developed
within this scheme the effective land use may increase with up to 60%. Thе concept hydro-
ecological monitoring of the Middle Ob River floodplain has been developed on a platform
for the organization of scientifically based, regionally adapted, and ecologically regulated
nature management.
Author to express one's thanks of professor Tomsky State University, PhD A.M. Adam for
data presentation and cause in hard expedites.
4. References
Adam A. M. & Bolotnov, V. P., (1982). Analysis of Influence of the Spring High Water over the
Structure of the Bird Population in the Floodplain of the Middle Ob for the Purpose of
Nature Protection. Deposited in VINITI No. 1040-82 [in Russian].
Adam A. M., Bolotnov V. P., and Sekisova S. E., (2001), inProblems of Geography of Siberia
.TGU, Tomsk, , Issue 24, pp. 211-218 [in Russian].
Adam A. M. and Mamin R. G., Natural Resources and Ecological Safety of West Siberia
(POLTEKS, Moscow, 2000) [in Russian].
Adam A. M., Novoselov A. L., and Chenurnykh N. V. (2000), Ecological Problems of the
Regions of Russia (VINITI, Moscow,) [in Russian
Bolotnov V. P. , Sekisova,S. E., and Adam A. M., (2001). "Environment of Siberia, the Far
East, and the Arctic," in Selected Paper Presented at the International Conference ESFA
2001, Tomsk, Russia (International Research Center of Environmental Physics and
Ecology, Russian Academy of Science, 2001), pp. 348-361.
Gul'tyaev A. K., (1999) MATLAB 5.2.1. Imitation Modeling inWindows: Practical Manual
(KORONA print, Sankt-Petersburg) [in Russian].
D'yakonov V. P. ,Abramenkova I. V., and Kruglov V. V., (2001). MATLAB 5.2.1 with Bump
Packs (Knowledge, Moscow) [in Russian].
Forrester J., (1978).World Dynamics (Nauka, Moscow) [Russian translation].
Ravkin, Yu. S., (1984). Spatial Organization of the Bird Population in the Forest Zone (West and
Central Siberia) (Nauka, Siberian Branch, Novosibirsk,) [in Russian].
Shepeleva L. F. (1986), Ekologiya, No. 2, 3.
Vartapetov D. G., (1984). Birds of the Taiga Interfluves of the West Siberia (Nauka, Siberian
Branch, Novosibirsk,) [in Russian].
Yudkin V. A., Ravkin Yu. S., Blinov V. N., et al., (1987). Spatial-Temporal Dynamics of the Fauna
(Birds and Small Mammals) (Nauka, Siberian Branch, Novosibirsk,) [in Russian].
5
1. Introduction
Constructing models, comparing their predictions with observations, and trying to
improve them, constitutes the core of the scientific approach to understanding complex
systems like large river basins (Even et al., 2007). These processes require manipulation of
huge historical data sets, which might be available in different formats and from various
stakeholders. The challenge is then to first pre-process the data to similar lengths, with
minimal loss of integrity, before manipulating it as per initial objectives. In the Upper and
Middle Vaal Water Management Areas (WMAs) of the Vaal River, bounded by Vaal dam
outlet and Bloemhof dam inlet, the overall objective of on-going research is to model
surface raw water quality variability in order to predict cost of treatment to potable water
standard. This paper reports on part of the overall research. Its objective was to show
how a huge and non-consistent water quality data set could be downsized to manageable
aspects with minimal loss of integrity. Within that scope, challenges were also
highlighted.
One of the more important forms of knowledge extraction is the identification of the more
relevant inputs. When identified, they may be treated as a reduced input for further
manipulation. In water quality data analysis, data collection, cleaning and pre-processing
are often the most time-consuming phases. All inputs and targets have to be transferred
directly from instrumentation or from other media, tagged and arranged in a matrix of
vectors with the same lengths (Alfassi et al., 2005). If vectors have outliers and/or missing
values these have to be identified for correction or to be discarded. More complex
mathematical correlations are sometimes employed to identify redundant, co-linear inputs,
or inputs with little information content (Alfassi et al., 2005).
Sources and sinks of variables in hydrodynamics, also known as forcing functions, are the
cause of change in water quality (Martin et al., 1998). To capture intermediate scale
processes that are spotty in spatial extent, extensive sampling and averaging of the
calibration data over sufficient spatial scales is done to capture that condition over time.
Although many water constituents are non-conservative in nature, a few conservative ones
that approach ideal behaviour under limited conditions, could be used for modelling and
calibration.
72 Scientific and Engineering Applications Using MATLAB
The study area is a major focus of modelling and pollution tracing in the Vaal basin, South
Africa, (Dzwairo et al., 2010b, Cloot and Roux, 1997, DWAF, 2007, Gouws and Coetzee,
1997, Naicker et al., 2003, Pieterse et al., 1987, Stevn and Toerien, 1976, Dzwairo et al., 2010a,
Dzwairo and Otieno, 2010, Herold et al., 2006).
Data sets spanning many years have been collected by various stakeholders including the
Department of Water Affairs (DWA) and Water Boards which treat bulk water for potable
use. For management of the basin as a whole these data sets come handy but the major
challenge is collating them into uniform and useable data, while noting that the different
stakeholders monitor selected parts of the basin for their own specific purposes. Some
sampling points might be dropped off or new points picked up as emerging pollution
threats require tracing and monitoring in order to mitigate effects. Still a useable data set
has to be constructed to monitor pollution and other threats, in addition to informing and
alerting decision makers regarding environmental and human health issues. This paper
shows how inconsistent and scattered data sets from 13 monitoring points were pre-treated
and downsized to SO42- inter-relationships. SO42- is a very important parameter in surface
water quality variability in this region because of the existence of gold and coal mining
activities. Threats from acid mine drainage are real.
2. Study area
The study area as indicated in Fig. 1 shows spatial relationships of the sampling points
located on VR and its tributaries as follows: B1-B10 on Blesbokspruit River (BR); K10-K10,
K6-K25 and K9-K19 on Klip River (KR); K12-N8 on Natalspruit River (NR); K1-R2 on
Withokspruit River, which is a tributary of Rietspruit River (RR); K3-R3 on another tributary
of RR; K2-R1 and K4-R4 on RR; S1-S1 and S4-S2 on Suikerbosrant River (SR); and V7-VRB37
and V9-VRB24 on Vaal River (VR).
B1-B10 K10-K10 K12-N8 K1-R2 K2-R1 K3-R3 K4-R4 K6-K25 K9-K19 S1-S1 S4-S2 V7-VRB37 V9-VRB24
1 2 3 4 5 6 7 8 9 10 11 12 13
Fig. 1. Monitoring points in study area bounded, by the two dams.
% Load the data with lots of missing dates. Note that in this example
% missing dates are not represented by NaN but are left out completely
>>[data,textdata] = xlsread('book.xls');
% Convert the text date to date numbers (you may have to change the date
% format depending on how your dates appear in Excel)
>>dates = datenum(textdata,'mm/dd/yyyy');
>>plot(dates,data,'LineStyle','none','Marker','o')
>>datetick('x')
% Create a new date series starting at the first date in dates and
% ending at the last but with every date in-between
>>newDates = dates(1):dates(end);
>>newData = interp1(dates,data,newDates,'cubic');
>>stringDates = cellstr(datestr(newDates));
>>newDates = dates(1):dates(end);
%Run the tic toc (3 instructions below at once by copying and pasting, it should
give elapsed time as eg 0.305720 seconds)
>>tic
newColumnData = interp1(dates,columnData,newDates,'cubic');
toc
%In a new figure, plot both the new data and the existing data
figure
>>plot(newDates,newColumnData,dates,columnData,'LineStyle','none','Marker','o')
>>datetick('x')
>tic
newColumnData = interp1(dates,columnData,newDates,'cubic');
toc
Warning: NaN found in Y, interpolation at undefined values will result in undefined values.
In interp1 at 178
Warning: All data points with NaN in their value will be ignored.
In polyfun\private\chckxy at 103
In pchip at 59
In interp1 at 283
Table 5. NaN.
Downsizing Water Quality Data for River Basin Management –
Focussing on Sulphate: Vaal River, South Africa Case Study 77
Another common error was that of a misplaced decimal point or full stop during data
capture (Table 6). Matlab would not be able to manipulate this entry for interpolation
because it was not a value. A duplicated or non-formatted date would also present an error
that would require debugging before a complete interpolated data set could be obtained.
These, among other similar errors, required manual debugging through a whole data set,
each a 2526 x28 matrix. With a perfect matrix, an interpolation took a fraction of a second.
4. Results
Initial inspection indicated that the data exhibited gross temporal inconsistency. Sampling
dates did not match, in addition to missing values. Table 7 shows the interpolated data for
points Z and Y for 5 to 21 July 2004.
A full length raw data set for Z (2003 to 2009), shown in Fig. 2, was interpolated and
graphed in Fig. 3, for only 4 out of the 27 variables, that is, Chl-α, COD, EC and DOC, to
reduce congestion and enhance clarity to the cubic interpolation concept.
Whereas Fig. 2 showed a legend with 4 data sets, Fig. 3’s legend included the interpolated
data, colour-coded for clarity. IChla, Icod, Iec and Idoc (IChl-α, ICOD, IEC and IDOC)
represented the interpolations of the 4 variables used. Daily interpolation was chosen for this
study because after interpolation, any other data interval, for example monthly or yearly
variation, could be computed without repeating the time-consuming interpolation process.
| cn_ ec do fc Hg Cl_ f_
-------------+---------------------------------------------------------------
cn_ | 1.0000
ec | 0.0908* 1.0000
do | -0.0106 0.0112* 1.0000
fc | 0.0014 0.0217* 0.0141* 1.0000
Hg | -0.0523* -0.1087* 0.0110 -0.0594* 1.0000
Cl_ | 0.0783* 0.8699* 0.0039 0.0062 -0.0192* 1.0000
f_ | -0.0053 0.1819* -0.0404* 0.0239* -0.1666* 0.0259* 1.0000
no2_ | -0.0708* -0.1365* 0.1629* 0.0809* 0.1839* -0.0458* -0.0787*
no3_ | -0.0628* 0.1223* 0.1033* 0.0658* 0.1916* 0.0876* 0.0115*
so42_ | 0.0961* 0.8720* -0.0064 0.0288* -0.2013* 0.7273* 0.2798*
Low_Hg | -0.0009 0.2998* 0.0450* -0.0260* -0.2516* 0.1762* 0.3496*
Mn | 0.0147* 0.3936* -0.0102 0.0668* -0.1783* 0.1815* 0.2316*
pH | 0.0290* -0.4242* 0.0481* -0.0856* 0.1456* -0.1382* -0.3480*
po43_ | -0.0367* -0.0858* 0.0283* 0.0418* 0.1250* -0.0193* -0.0683*
s | 0.0807* 0.8861* -0.0176* 0.0226* -0.1974* 0.7435* 0.2593*
ss | -0.0302* -0.2024* -0.0336* 0.0138 0.0350* -0.1852* -0.0387*
Temp | -0.0120* -0.0369* -0.0424* 0.0201* -0.0948* -0.0544* 0.0481*
T_Silica | -0.0343* 0.1377* -0.0693* 0.0422* -0.1797* -0.0889* 0.2674*
Turb | -0.0434* -0.2525* -0.0862* 0.0284* -0.0893* -0.2899* 0.0213*
nh4_ | 0.0267* 0.3493* -0.0444* 0.2118* -0.0952* 0.2378* 0.1670*
Chla | 0.0039 0.0918* 0.1341* -0.0320* 0.0218 0.1432* 0.0204*
cod | -0.0546* -0.2345* -0.0950* 0.0367* -0.2205* -0.1833* -0.1091*
doc | -0.0661* -0.4022* -0.0080 -0.0702* 0.0607* -0.2446* -0.1826*
Mo | -0.0172* -0.0089 0.0123* 0.0099 -0.0743* 0.0042 0.1316*
Si | -0.0335* 0.1380* -0.0697* 0.0420* -0.1789* -0.0880* 0.2640*
p | -0.0621* -0.1345* 0.0126* 0.0885* 0.1870* -0.0679* -0.0701*
80 Scientific and Engineering Applications Using MATLAB
| cod doc Mo Si p Fe
-------------+------------------------------------------------------
cod | 1.0000
doc | 0.5436* 1.0000
Mo | 0.0334* 0.0810* 1.0000
Si | -0.0168* -0.2441* -0.0451* 1.0000
p | 0.0381* 0.1008* 0.0430* 0.0570* 1.0000
Fe | -0.0369* -0.1302* -0.0176* 0.3519* -0.0767* 1.0000
--------------------------------------------------------------------------
Factor | Eigenvalue Difference Proportion Cumulative
----------------+------------------------------------------------------------
Factor1 | 5.82041 3.19894 0.5510 0.5510
Factor2 | 2.62148 0.50078 0.2482 0.7992
Factor3 | 2.12070 1.29933 0.2008 1.0000
-----------------------------------------------------------
Variable | Factor1 Factor2 Factor3 | Uniqueness
-------------+------------------------------+--------------
cn_ | | 0.9977
ec | 0.6603 | 0.4260
do | | 0.9881
fc | | 0.9666
Hg | -0.4816 | 0.7544
Cl_ | 0.7176 | 0.1997
f_ | | 0.9921
no2_ | 0.5019 | 0.7768
no3_ | 0.8243 | 0.3693
so42_ | 0.8206 | 0.2361
Low_Hg | 0.6888 | 0.6217
Mn | 0.7274 | 0.5483
pH | -0.4832 | 0.6090
po43_ | | 0.9908
s | 0.8318 | 0.2598
ss | 0.8475 | 0.3456
Temp | 0.3315 | 0.8679
T_Silica | 0.6666 | 0.2333
Turb | 0.8739 | 0.2462
nh4_ | 0.7095 | 0.5037
Chla | | 0.8587
cod | 0.6745 0.4000 | 0.5787
doc | 0.7211 0.3964 | 0.4579
Mo | 0.4133 | 0.8677
Si | 0.6684 | 0.2326
p | | 0.9023
Fe | 0.6249 | 0.6065
-----------------------------------------------------------
(blanks represent abs(loading)<.33)
Table 10. Rotated factor loadings (pattern matrix) and unique variances.
EC and Cl-, together with FC, Hg, F-, NO3-, Low_Hg, Mn, pH, S, SS, Temp, T_Silica, Turb,
NH4+, COD, Si, P and Fe, were good predictors for SO42- concentration, and the fitted model
explains 82% of the total variation (Table 12).
------------------------------------------------------
Variable | Factor1 Factor2 Factor3 | Uniqueness
------------------------------------------------------
ec | 0.6603 | 0.4260
Hg | -0.4816 | 0.7544
Cl_ | 0.7176 | 0.1997
no2_ | 0.5019 | 0.7768
no3_ | 0.8243 | 0.3693
so42_ | 0.8206 | 0.2361
Low_Hg | 0.6888 | 0.6217
Mn | 0.7274 | 0.5483
pH | -0.4832 | 0.6090
s | 0.8318 | 0.2598
ss | 0.8475 | 0.3456
Temp | 0.3315 | 0.8679
T_Silica | 0.6666 | 0.2333
Turb | 0.8739 | 0.2462
nh4_ | 0.7095 | 0.5037
cod | 0.6745 | 0.5787
doc | 0.3964 | 0.4579
Mo | 0.4133 | 0.8677
Si | 0.6684 | 0.2326
Fe | 0.6249 | 0.6065
-----------------------------------------------------------
(blanks represent abs(loading)<.33)
------------------------------------------------------------------------------
so42_ | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
cn_ | -22.32404 18.52691 -1.20 0.228 -58.64195 13.99386
ec | .3736444 .0227941 16.39 0.000 .3289616 .4183271
do | .0131522 .0926716 0.14 0.887 -.1685098 .1948143
fc | .0000566 .0000189 2.99 0.003 .0000195 .0000938
Hg | -89.09861 11.70687 -7.61 0.000 -112.0473 -66.14989
Cl_ | .7573463 .042237 17.93 0.000 .67455 .8401425
f_ | 32.3612 8.280861 3.91 0.000 16.12841 48.59399
no2_ | -10.90126 13.10631 -0.83 0.406 -36.59327 14.79075
no3_ | 3.180277 1.003154 3.17 0.002 1.213816 5.146738
Low_Hg | -4.527516 .8473181 -5.34 0.000 -6.188495 -2.866536
Mn | 51.43273 4.405735 11.67 0.000 42.79626 60.06919
pH | -7.478322 2.569807 -2.91 0.004 -12.51586 -2.440786
po43_ | .8106866 .7992836 1.01 0.310 -.7561315 2.377505
s | 1.743953 .0246683 70.70 0.000 1.695596 1.79231
Downsizing Water Quality Data for River Basin Management –
Focussing on Sulphate: Vaal River, South Africa Case Study 83
| Summary of so42_
Sample_ID | Mean Std. Dev. Freq.
------------+------------------------------------
B1-B10 | 405.26118 140.67122 2526
K1-R2 | 66.18701 115.52301 2526
K10-K10 | 120.27818 58.483346 2526
K12-N8 | 303.80768 116.03529 2526
K2-R1 | 1128.8242 815.12126 2526
K3-R3 | 121.64965 170.8744 2526
K4-R4 | 1123.08 607.58752 2526
K6-K25 | 172.05588 44.633777 2526
K9-K19 | 163.85514 45.159634 2526
S1-S1 | 21.228942 11.581847 2526
S4-S2 | 346.77498 144.27252 2526
V7-VRB37 | 159.3354 44.584895 2526
V9-VRB24 | 154.30907 45.776534 2526
------------+------------------------------------
Total | 329.7421 462.44325 32838
Analysis of Variance
Source SS df MS F Prob > F
------------------------------------------------------------------------
Between groups 4.1391e+09 12 344925487 3926.94 0.0000
Within groups 2.8832e+09 32825 87835.795
------------------------------------------------------------------------
Total 7.0223e+09 32837 213853.757
(Sidak)
Row Mean-|
Col Mean | B1-B10 K1-R2 K10-K10 K12-N8 K2-R1 K3-R3
---------+------------------------------------------------------------------
K1-R2 | -339.074
| 0.000
K10-K10 | -284.983 54.0912
| 0.000 0.000
K12-N8 | -101.453 237.621 183.529
| 0.000 0.000 0.000
K2-R1 | 723.563 1062.64 1008.55 825.017
| 0.000 0.000 0.000 0.000
K3-R3 | -283.612 55.4626 1.37148 -182.158 -1007.17
| 0.000 0.000 1.000 0.000 0.000
K4-R4 | 717.819 1056.89 1002.8 819.272 -5.7442 1001.43
| 0.000 0.000 0.000 0.000 1.000 0.000
K6-K25 | -233.205 105.869 51.7777 -131.752 -956.768 50.4062
| 0.000 0.000 0.000 0.000 0.000 0.000
K9-K19 | -241.406 97.6681 43.577 -139.953 -964.969 42.2055
| 0.000 0.000 0.000 0.000 0.000 0.000
S1-S1 | -384.032 -44.9581 -99.0492 -282.579 -1107.6 -100.421
| 0.000 0.000 0.000 0.000 0.000 0.000
S4-S2 | -58.4862 280.588 226.497 42.9673 -782.049 225.125
| 0.000 0.000 0.000 0.000 0.000 0.000
V7-VRB37 | -245.926 93.1484 39.0572 -144.472 -969.489 37.6857
| 0.000 0.000 0.000 0.000 0.000 0.000
V9-VRB24 | -250.952 88.1221 34.0309 -149.499 -974.515 32.6594
| 0.000 0.000 0.004 0.000 0.000 0.007
Row Mean-|
Col Mean | K4-R4 K6-K25 K9-K19 S1-S1 S4-S2 V7-VRB37
---------+------------------------------------------------------------------
K6-K25 | -951.024
| 0.000
K9-K19 | -959.225 -8.20074
| 0.000 1.000
S1-S1 | -1101.85 -150.827 -142.626
| 0.000 0.000 0.000
S4-S2 | -776.305 174.719 182.92 325.546
| 0.000 0.000 0.000 0.000
V7-VRB37 | -963.745 -12.7205 -4.51974 138.106 -187.44
| 0.000 1.000 1.000 0.000 0.000
V9-VRB24 | -968.771 -17.7468 -9.54607 133.08 -192.466 -5.02633
| 0.000 0.929 1.000 0.000 0.000 1.000
Multivariate linear regression indicated that out of the 26 variables that could predict SO42- ,
only 20 were significant, accounting for 82% of the total variation of SO42-.
While correlation and regression provided linear relationships, factor analysis, on the other
hand, could be used for data reduction. Even though sometimes it is difficult to find a
common name to assign to a factor, still, based on these statistical approaches, individual
factors or elements within a factor could be further analysed as necessary, with minimal loss
of data integrity.
From one-way ANOVA, SO42- mean concentration values indicated that monitoring point
K2-R1 (1128.82±815 mg/L) was within the vicinity of the source of SO42-. Attenuation of the
variable was noted as its mean value decreased along the Rietspruit River at K4-R4 and then
Klip River at K6-K25 and K9-K19, before Klip River discharged into the Vaal River. From
monitoring point B1-B10 (also close to a source of SO42-), another established route was
through S4-S2, before Suikerbosrant River discharged into the Vaal River upstream of the
Klip River. Surface raw water containing high levels of SO42- was not draining via K1-R2
and S1-S. Based on SO42- mean concentration values only and for management purposes,
K1-R2 and S1-S could be left out of the monitoring programme, saving on financial
resources. Comparison of SO42- by sample_ID showed that K6-K25, K9-K19, V7-VRB37 and
V9-VRB24; K10-K10 and K3-R3; and K2-R1 and K4-R4, were significantly similar.
The major challenge was pre-processing of the non-consistent water quality data over the 7
years. Non-consistent data was as a result of missing data, largely where some of the
stakeholders dropped or established some water quality variables and monitoring points
over the years as monitoring prioritizations changed because of new and emerging
pollution threats. The challenge of insufficient and inconsistent data for water quality
modelling remains a limitation in the formulation of good and practically useable models.
However, interpolations and correlations, including factor analysis and regression, could
help build better data sets, especially for pollution trending in river basin management. This
could be used to support large-scale public decisions.
6. Acknowledgement
The financial assistance of the South African Department of Science Technology (DST) is
hereby acknowledged. Opinions expressed and conclusions arrived at, are those of the
authors and are not necessarily to be attributed to the DST. The authors would also like to
thank Tshwane University of Technology for hosting and co-funding this research. DWA,
the Water Research Commission, Rand Water Board (co-funding), Midvaal Water Company
and Sedibeng Water, are also sincerely acknowledged, especially for providing very
valuable and vital data.
7. References
Alfassi, Z. B., Boger, Z. & Ronen, Y. (2005). Statistical treatment of analytical data, Oxford,
Blackwell Science Ltd, 0-632-05367-4, CRC Press, Australia.
Cloot, A. & Roux, G. L. (1997). Modelling algal blooms in the middle Vaal River: a site
specific approach. Water Research, 31, 2, 271-279, 0043-1354.
DWAF (2007). Integrated water quality management plan for the Vaal River system.
Pretoria, South Africa.
86 Scientific and Engineering Applications Using MATLAB
Dzwairo, B. & Otieno, F. A. O. (2010). Integrating quality and cost of surface raw water:
Upper and Middle Vaal Water Management Areas South Africa. Water Science and
Technology: Water Supply 10, 2, 201–207, 1606-9749.
Dzwairo, B., Otieno, F. A. O. & Ochieng', G. M. (2010a). Making a case for systems thinking
approach to integrated water resources management (IWRM). International Journal
of Water Resources and Environmental Engineering, 1, 5, 107-113 2141-6613.
Dzwairo, B., Otieno, F. A. O., Ochieng', G. M. & Letsoalo, M. A. (2010b). Downsizing water
quality data for river basin management – Focussing on Sulphate: Vaal River, South
Africa. Proceedings of the 11th WaterNet/WARFSA/GWP-SA Symposium: ‘IWRM for
National and Regional Integration: Where Science, Policy and Practice Meet. Elephant
Hills Hotel, Victoria Falls, Zimbabwe. 27 October - 29 October, 2010.
Even, S., Billen, G., Bacq, N., Théry, S., Ruelland, D., Garnier, J., Cugier, P., Poulin, M., Blanc,
S., Lamy, F. & Paffoni, C. (2007). New tools for modelling water quality of
hydrosystems: An application in the Seine River basin in the frame of the Water
Framework Directive. Science of The Total Environment, 375, 1-3, 274-291.
Gouws, K. & Coetzee, P. P. (1997). Determination and partitioning of heavy metals in
sediments of the Vaal Dam System by sequential extraction. Water SA, 23, 3, 217-
226, 0378-4738.
Herold, C. E., Le Roux, P. J., Nyabeze, W. R. & Gerber, A. (2006). WQ2000 Salinity Model:
enhancement, technology transfer and implementation of user support for the Vaal
system. Umfula Wempilo Consulting. Pretoria, South Africa.
Ho, R. (2006). Handbook of univariate and multivariate data analysis and interpretation with SPSS.,
Florida, Chapman and Hall/CRC: Tailor and Francis Group, 1584886021.
Martin, J. L., Mccutcheon, S. C. & Martin, M. L. (1998). Hydrodynamics and Transport for Water
Quality Modeling Taylor & Francis, Inc, 978-0873716123.
Naicker, K., Cukrowska, E. & Mccarthy, T. S. (2003). Acid mine drainage arising from gold
mining activity in Johannesburg, South Africa and environs. Environmental
Pollution, 122, 1, 29-40, 0269-7491.
Ochse, E. (2007). Seasonal rainfall influences on main pollutants in the Vaal River barrage reservoir:
a temporal-spatial perspective. Magister Artium MA, University of Johannesburg.
Pieterse, A., Roos, J., Roos, K. & Pienaar, C. (1987). Preliminary observations on cross-
channel and vertical heterogeneity in environmental and algological parameters in
the Vaal River at Balkfontein, South Africa. Water SA, 12, 4, 173-184, 0378-4738 .
Stevn, D. J. & Toerien, D. F. (1976). Eutrophication levels of some South African
impoundments. IV. Vaal dam. Water SA, 2, 2, 53-57.
6
1. Introduction
There is a growing concern on the capacity of water utilities via governmental regulatory
agencies regarding potential optimization and reliability for water distribution network.
Generally, water distribution networks comprise about 60% of the total budget for a
complete framework of a water supply system. According to this fact, achieving an
optimum solution for water distribution networks as models’ outcome of reliability-based
optimization design has become the great concern to save considerable amount of allocated
budget. During the last decade many authors were interested in studying optimization and
reliability for water distribution networks that include solving non-linear hard problem of
the network hydraulic equations. The optimization and reliability models of water
distribution networks have number of varieties in studying aspects that include efficiency,
accuracy, different sizes/scales of networks, and the consumed run time to define the
optimum solution. During the current decade, considerable amount of attention has been
given to reliability of water distribution networks in conjunction with the optimization to
achieve maximum benefits with the minimum cost. This concern has been extended to cover
the risk management for water distribution networks as a way to embark on facing the
shortage of water resources all over the world or improving asset management programs.
The main objective of this chapter is to develop standalone model divided into four sub-
models using MATLAB environment programming language. The developed model and its
corresponding sub-models would acquaint an optimum solution for a given water
distribution network that achieve both least cost design and reliability based optimization
design in the mean time. The main model is called RELOPT and can be used as a tool to
implement: modeling reliability-based optimization design, deterioration analysis of water
pipe networks, risk analysis and assessment, and decision support system. RELOPT is
integrated with four sub-models those are: optimization search engine model that is based
on a new technique driven from Genetic Algorithms approach is called Linear Adaptive
Genetic Algorithm (LAGA); pre-estimation optimization model that is based on Average
Gradient Method (AGM) to accelerate the process of the optimization search engine;
reliability model that is based on load resistance concept for calculating system reliability.
Through this chapter the number of subjects will be discussed those are: background of
water distribution systems, definition of problem in statement, main objectives, history of
optimization and techniques, history of reliability and techniques, proposed optimization
technique, the advantage of the new optimization technique, proposed reliability evaluation
88 Scientific and Engineering Applications Using MATLAB
2. Chapter objectives
The objectives of this chapter have one common target that is to define the reliability-based
optimization design for a water distribution network using modelling technique of
MATLAB programming language. The following tasks have to be achieved:
1. Acquainting optimum least-cost design for water distribution networks using new
efficient and time consumed method.
2. Define risk components for water distribution networks.
3. Define the most critical components of water distribution networks that affect the level
of serviceability under different cases of operation (i.e. define level of service under
risk).
4. Analysis, evaluation and treating reliability for water distribution networks.
5. Define the reliability of water distribution network over a given period of time.
6. Define the optimum solution of water distribution network that achieve the optimum
lease-cost design and certain accepted reliability in one time (reliability-based
optimization).
7. Develop stand alone reliability-based optimization model comprising all the above
mentioned objectives
3. Previous studies
Solving the hydraulic equations for water distribution networks is a constrained non-linear
hard programming problem (CNLHP) due to the nature of the non-linearity of the decision
variables such as pipe diameters. For a given water distribution network, huge number of
solutions could be selected through a range of the decision variables to select the best
solution which arise the problem to be combinatorial optimization problems (Gupta and
Kapoor 1994). Hamdy A. (1997) stated that some mathematical models may be so complex
that is impossible to solve them by any of the available optimization algorithm and such
cases heuristics are used instead of mathematical models to search for a good solution near
the best or the optimum one. The advantage of heuristics over an exact optimization
algorithm is that it is usually much faster to execute. Dorigo and Thomas (2004) stated that
recently, many researchers have focused their attention on a new class of algorithms called
meta-heuristics. A meta-heuristic is a set of algorithmic concepts that can be used to define
heuristic methods applicable to a wide set of different problems. The use of meta-heuristics
has significantly increased the ability of finding very high quality solutions to hard,
practically relevant combinatorial optimization problems in reasonable time. Intelligent
algorithms models are becoming essential tool to solve such non-linear hard problems for
water distribution networks. Number of intelligent algorithms had been developed based
on the meta-heuristic concept such as Simulating Annealing (SA), Tabu Search (TS), Guided
Local Search (GLS), Greedy Randomized Adaptive Search Procedure (GRASP), Iterated
Local Search (ILS), Evolutionary Computation (EC), Scatter Search, and Ant Colony
Optimization (ACO).
Modelling Reliability Based Optimization Design for Water Distribution Networks 89
4. Model components
Getting the reliability-based optimization design for water distribution networks requires
searching among a number of available population set of solutions, thus; RELOPT model
consists of the following components (Moneim AM, 2009):
1. Hydraulic solver EPANET: consists of the dynamic libraries that are required to be
called by MATLAB program for hydraulic analysis.
2. Pre-estimation model (AGM): this sub-model provides the lower and upper bounds
that are required for OPTWNET to start optimization search process.
3. Optimization model (OPTWNET): defines the optimum solution using LAGA.
4. LAGA automatic search engine module.
5. Reliability model (RELWNET): this model is connected with three sub-models those
are: minimum cut-sets model; Generic Expectation Function model; and reliability
calculation model. The model passes the final calculated reliability to the main model
RELOPT.
used by LAGA is limited to 5 and the population size is limited to 10. Experiments have
shown that these limits are enough for LAGA to define the global optimization solution for a
water network. OPTWNET passes the resulting optimized decision variables to the main
model which passes by its role to the reliability model RELWNET.
5. Optimization technique
Heuristics from Nature approach has been presented by (Colorni et al. 1992a, 1992b)
presented methods that depend on as a non-derivative optimization method. The authors
stated that Heuristics derived from Nature algorithms considered as a border between
Operation Research (OR) and Artificial Intelligence (AI). These algorithms take inspiration
from physics, biology, social, and use a certain amount of repeated trials for the Non-
Programming Hard Combinatorial Optimization Problem (NPHCOP). Heuristics are
obtained by one of the following application methods:
1. Using certain amount of repeated trials;
2. Employing one or more "agents" such as neurons, particles, chromosomes, ants, genetic
algorithms, and so on;
3. Operating (in case of multiple agents) with a mechanism of competition-cooperation;
4. Embedding procedures of self-modifications of the heuristics parameters or of the
problem representation.
0.5
Pc = (g 1)+1 (1)
M 1
Where: g= number of the current generation,
M= total number of generation.
Unlike the crossover probability is the mutation probability which is not needed at the
beginning of the optimization search as the population is very distinct. As the generation is
progressing the solution starts to come slightly closer to the optimum, the mutation
probability is then come to the view picture and starts from 0.005 and increased linearly up
to 0.5.
0.005
Pm = (g 1) (2)
M 1
range and substitute directly in the objective function to evaluate the corresponding selected
input variables. LAGA has been adopted to select the population for network solution
(string of pipe diameters) and linked to hydraulic solver EPANET to evaluate each network
solution according to predefined cost function. The second problem is that the search space
for given pipe diameter limit was very large and that slowed down the convergence process.
This problem had been eliminated by applying pre-estimation for the upper and lower
bounds of diameter population as being explained in the forthcoming section 3.14. Figure
3.5 illustrates the flow chart for LAGA optimization technique. As being introduced in this
chapter, the problem of water distribution networks need evaluation of each selected
network solution by linking to hydraulic solver such as EPANET. This process is affecting
the consumed run time for specific problem and as a result the expected run time to get the
final solution will be much longer. The following advantages of LAGA technique are
helping to reduce and accelerate the run time:
1. LAGA uses dynamic concept for the main parameters of GA process that allow
changing the values of both cross over and mutation probability at the start of each
generation. This mechanism is accelerating the searching process to allocate an
optimum solution in shorter time than the traditional Genetic Algorithm;
2. LAGA reduces the population size as the optimization process tends to catch on the
optimum solution. Conversely, LAGA increases the population size as the optimization
process could not catch up an optimized solution. This mechanism adds another time
consuming facility in case an optimum solution has been identified.
intended to a value greater than zero. The objective function for water distribution networks
is usually used as the cost function comprising the pipe lengths and diameters. The cost
function can take the following general form:
n
MIN(F)= C i *D i *l i +Fp +Fv (4)
i=1
Where F = total cost of the current network and MIN is referring to minimization.
Fp =penalty cost due to pressure violation.
Fv = penalty cost due to velocity violation
C= unit cost for each pipe diameter category
D= pipe diameter
l= pipe length
i=pipe number
n=denotes to total pipe number within the network
F
Fp = (R p -Pmin ) (5)
Rp
If the minimum allocated pressure of the current solution is greater than the desirable
minimum pressure the resulting penalty cost will have negative value which will decrease
the total corresponding cost function. Conversely, if the minim allocated pressure of the
current solution is less than the desirable minimum pressure, the resulting penalty cost will
have positive value which will increase the corresponding cost function.
1. Case of minimum allocated pressure Pmin is less than minimum desirable pressure R p
and the difference does not exceed 10%:
F
Fp = (R p -Pmin ) *1000 (6)
Rp
2. Case of more than 10% difference between the actual minimum pressure and the
minimum desirable pressure and the solution is invisible:
F
Fp = |R p -Pmin|*1000 (7)
Rp
3. Case number 2 and 3 will eliminate the solutions that violate about the optimum
solution even if it is a feasible solution. The factor 1000 (this factor can be increased to
any other value) will increase the cost value to an imaginary amount which secures
disappearing of any violated solutions that exceed 10 % of the desirable minimum
pressure.
96 Scientific and Engineering Applications Using MATLAB
-1*F
Fv = |Va -V|
r (8)
Vr
Where Vr = desirable maximum velocity.
Va = Actual maximum velocity of the current solution.
If the actual maximum velocity of the current solution is greater or less than the desirable
maximum velocity within a limit of 10% of the resulting penalty cost will have negative
value which will decrease the total corresponding cost function.
2. Case of more than 10% difference between the actual maximum velocity and the
maximum desirable velocity and the solution is visible:
-1*F
Fp = |Vr -Va|*1000 (9)
Vr
Case number 2 will eliminate the solutions that violate about the optimum solution even if it
is a visible solution. High velocities might cause pipe corrosion and high power loss. The
factor (1000) will increase the cost value to an imaginary amount which securing
disappearing of any violation solutions that exceed 10 % of the desirable maximum velocity.
essentially needed to develop a reliability model. The system Reliability RS can be calculated
after the probability of system failure can be calculated. MATLAB built-in probability
functions can be used to generate any number of historical data for both the load that is
carried by a pipe and the resistance that is can be calculated using Hazen-William equation.
S
PF= P(C k ) (11)
k=1
Where PF = Probability of failure of the whole system
P(C k ) = Probability of failure for the minimum cut set number k
S = total number of the minimum cut set.
It is easier to define the failure components and calculate the failure probability of each one,
then the overall system failure cab be calculated simply according to the following sections.
M
RS=1-P=1- (1-HA)*PF (13)
i
Where: RS = Overall system reliability
P = Overall system failure
98 Scientific and Engineering Applications Using MATLAB
8. Application example
A simple 2-loop network represented in Figure 4 consists of 7 nodes, one reservoir and 8
pipes have been obtained from the literature (Alperovits & Shamir, 1977). All pipes are
equally length 1000 m and the Hazen-Williams coefficient of friction is assumed to be 130
for all pipes. The minimum required residual pressure for all demand nodes is 30m. The
unit cost for pipe diameters is represented in Table 1 while Table 2 provides the data of
nodes’ demands and elevations.
Diameter
1 2 3 4 6 8 10 12 14 16 18 20 22 24
(In.)
Unit
2 5 8 11 16 23 32 50 60 90 130 170 300 550
Cost
Table 1. Unit Cost per meter length of pipe.
Node
1 2 3 4 5 6 7
Number
Demand
-1120 100 100 120 270 330 200
(m3/hr)
Elevation
210 150 160 155 150 165 160
(m)
Table 2. Node number, demands and elevation.
45
40
35
Residual Pressure (m)
30
25
20
15
10
0
0.95 0.96 0.97 0.98 0.99 1
System Reliability
Fig. 5. Reliability versus Residual Pressure for 2-loop Network.
6
x 10
4.5
3.5
2.5
Cost
1.5
0.5
0
0.95 0.96 0.97 0.98 0.99 1
System Reliability
N
C= 6.85*e0.1312D *L (14)
i=1
Where: D i = Pipe diameter in inches for pipe and L i is the pipe length in meters.
The results indicated that relationship between the cost and the reliability is directly
proportional which is logically expected. The point of issue is that a characteristic equation
can express the relationship between the cost and reliability and accordingly the cost can be
calculated for a given value of reliability. This can help when specific allocated budget needs
to be adopted for a given water network. Fig. 7 illustrates the relationship between the cost
and system reliability for the 2 loop water network.
102 Scientific and Engineering Applications Using MATLAB
9. Conclusion
In this chapter new model for reliability-based optimization has been developed using
MATLAB programming language to define an optimum and reliable design for water
distribution networks. The results indicated that the model RELOPT provides acceptable
level of confidence when has applied to reliability-based optimization problems. RELOPT
consist of four models those are: Pre-estimation mode called Average Gradient Model
(AGM); hydraulic solver EPANET; optimization model OPTWNET based on new genetic
algorithm concept called Linear Adaptive Genetic Algorithm (LAGA); and Reliability model
RELWNET based on load resistance concept for reliability evaluation. The minimum cut sets
method has been adopted to define the risk components for complex systems such as water
distribution networks. Generic Expectation Functions (GEF) has been adopted to define the
probability distribution for the difference between the load and resistance and hence the
failure probability has been calculated to calculate the system reliability. It has been
illustrated through that achieving reliability-based optimization represents the optimum
solution for water distribution networks when considered as constraint during the
optimization search process. Optimization of design-based cost is normally securing the
function of the water distribution networks during normal operation while reliability-based
optimization is normally securing the function of the water distribution network under risk.
Giving that the risk has been occurred is an important approach to apply during the design
of water distribution networks especially for large scale networks. Optimization cost for
water networks could save nations' budget while designing a new network but shall not
securing strategic decisions under certain risks or failures. On another point of view, while
Modelling Reliability Based Optimization Design for Water Distribution Networks 103
certain failure has occurred, budget may exceed normal levels to overcome such risks.
Reliability-based optimization helps decision makers to adopt pre optimum strategic
maintenance/operation programs or emergency scenarios to overcome lots of problems
caused by network failures. Increasing the current construction budget for optimization
based cost design of a water network to accommodate its reliability, will save more
investment in future giving that certain risks have occurred.
10. Acknowledgment
The author is so grateful to his professors Dr. Ahmed Moawad who has supported me along
my researches development with his value time and advices. The author also submit his
greetings to his professors Dr. Anas Al Molla, Dr. Ayman Al Salawy, Dr. Abdel Badee Salem
and Dr. Amgad El-Ansary who have given their value advices and encouraging through my
research and life of experience. Great thanks to everybody who has helped me with an
advice, time, word and/or technical data which with no doubt has a positive effect to
improve my work.
11. References
Alperovits E. and Shamir, U., (1977). Design of optimal water distribution systems. Water
Resources Research, Vol. 13(6), pp. 885-900.
Attia A. and Horacek P., (2001). Adaptation of genetic algorithms for optimization
problem solving. 7th Intern. Mendel conference on soft computing pp. 36-41, Brno,
Cizek Republic.
Colorni, A., M. Dorigo and V. Manniezzo (1992a). Distributed optimization by ant colonies.
In: Proceedings of the First European Conference on Artificial Life (ECAL-91)
(F.J. Varela and P. Bourgine, Ed.). The MIT Press. Cambridge MA. pp. 134–
142.
Colorni, A., M. Dorigo and V. Manniezzo (1992b). An investigation of some properties of an
ant algorithm". In: Parallel problem solving from nature, Vol 2. (R. M¨anner and B.
Manderick, Ed.). North-Holland. Amsterdam, pp. 509–520.
David A Coley, 1999. An Introduction to Genetic Algorithms for Scientists and Engineers.
World Scientific Publishing Company.
Dorigo, M., and Thomas, S., (2004). Ant colony optimization. MIT Press, Cambridge,
Massachusetts, London, England.
Gupta, I., Gupta, A., and Khanna, P. (1999). Genetic algorithm for optimization of water
distribution systems. Environmental Modelling & software Vol.-4, pp. 437-446.
Hamdy, A. T., 1997. Operation research an introduction. Sixth Edition, Prentice-Hall
International, Upper Saddle River, New Jersey.
Moneim A.M., Moawad A.K., Molla A., and Selawy A. (2010). RELWANET: Reliability
Evaluation Model for Water Distribution Networks. Australian Journal of Water
Resources, Volume 14 No. 1.
Su, Y. C., Mays, L. W., Duan, N. and Lansey, K. E. (1987). Reliability-based
optimizationmodel for water distribution system. Journal of Hydraulic
Engineering, ASCE, 114(12), 1539-1556.
104 Scientific and Engineering Applications Using MATLAB
1. Introduction
In cyber-physical systems (CPSs), embedded computing systems and communication capability
are used to streamline and fortify the operation of a physical system. Intelligent critical
infrastructure systems are among the most important CPSs and also prime examples of
pervasive computing systems, as they exploit computing to provide "anytime, anywhere"
transparent services. While the added intelligence offers the promise of increased utilization,
its impact must be assessed, as unrestricted cyber control can actually lower the reliability of
existing infrastructure systems.
As a practical example, water distribution networks (WDNs) are an emerging CPS domain.
Physical components, e.g., valves, pipes, and reservoirs, are coupled with the hardware
and software that support intelligent water allocation. An example is depicted in Fig. 1.
The primary goal of WDNs is to provide a dependable source of potable water to the
public. Information such as demand patterns, water quantity (flow and pressure head), and
water quality (contaminants and minerals) is critical in achieving this goal, and beneficial in
guiding maintenance efforts and identifying vulnerable areas requiring fortification and/or
monitoring. Sensors dispersed in the physical infrastructure collect this information, which is
fed to algorithms (often distributed) running on the cyber infrastructure. These algorithms
provide decision support to hardware controllers that are used to manage the allocation
(quantity) and chemical composition (quality) of the water. As WDNs become larger and
more complex, their reliability comes into question.
Modeling and simulation can be used to analyze CPS performability, as direct observation
of critical infrastructure is often infeasible. Accurate representation of a CPS encompasses
three aspects: computing, communication, and the physical infrastructure. Fundamental
differences exist between the attributes of cyber and physical components, significantly
complicating representation of their behavior with a single comprehensive model or
simulation tool. Specialized simulation tools exist for the engineering domains represented
in critical infrastructure, including power, water, and transportation. These tools have been
created with the objective of accurately reflecting the operation of the physical system, at high
spatial and temporal resolution. As is the case with specialized models of physical systems,
intelligent control is not reflected in these tools. Despite the existence of simulation tools for
cyber aspects such as computing and communication, differences in temporal resolution and
106
2 Scientific and Engineering Applications Using MATLAB
Will-be-set-by-IN-TECH
data representation and the lack of well-defined interfaces pose considerable challenges to
linking these simulation tools in a fashion that accurately represents the CPS as a whole.
In the first part of this chapter, we articulate the available simulation tools and the challenges
present in integrated simulation of CPS, where the goal is to accurately reflect the operation
and interaction of the cyber and physical networks that comprise the system. A solution is
presented for the CPS domain of intelligent WDNs. The proposed solution utilizes EPANET
to simulate the physical infrastructure of the water distribution network and Matlab to
simulate the cyberinfrastructure providing decision support. Communication between the
two simulators replicates the interactions between cyber and physical components of WDNs,
and facilitates the observation of physical manifestations of intelligent control decisions.
This communication between the simulators takes place without user intervention, as all
information relevant to each simulator has been identified and extracted from the output of
the other. Information flows from the physical simulator to the cyber simulator, replicating
the operation of sensors in the physical infrastructure. The cyber simulator processes this
data in Matlab, and provides decision support for water allocation, in the form of setting for
control elements in the physical infrastructure. This information is provided to the physical
simulator, which applies these settings. This process repeats for the duration of the simulation,
as it would in the actual operation of a CPS.
The second part of this chapter addresses computation in the CPSs, specifically, the role
of cyberinfrastructure in CPSs. We present an agent-based framework for intelligent
environmental decision support. Due to the flexibility of software agents as autonomous and
intelligent decision-making components, the agent-based computing paradigm is proposed
for surmounting the challenges posed by a) fundamental differences in the operation of
cyber and physical components, and b) significant interdependency among the cyber and
physical components. The environmental management domain used as a model problem
is water distribution, where the goal is allocation of water to different consuming entities,
subject to the constraints of the physical infrastructure. In the cyber-physical approach to this
problem, which is implemented by intelligent WDNs, the cyberinfrastructure uses data from
the physical infrastructure to provide decision support for water allocation. We adopt game
theory as the algorithmic technique used for agent-based decision support in an intelligent
Integrated Cyber-Physical
Integrated Cyber-Physical Simulation
Simulation of Intelligent of Intelligent
Water Distribution Networks Water Distribution Networks 1073
WDN. In this initial effort, our focus is on providing decision support for the quantity of water
allocated to each consuming entity. Game theory is a natural choice for complex resource
allocation problems such as water distribution, where hydraulic and physical constraints,
ethical concerns, and economic considerations should be represented. The investigation
of game theory as the computational algorithm for water quantity allocation is assisted
by Matlab, due to its powerful computational capability and ability to support advanced
techniques, such as distributed decision support algorithms. EPANET provides the data used
by the distributed computing algorithm to decide on water quantities.
In the third part of the chapter, we study the combination of game theory and the integrated
cyber-physical simulator, and investigate how different configuration of actuators based
on the game theory strategy can influence the malfunction of the purely physical WDN
in the EPANET. When the faults are injected into the physical infrastructure (represented
by EPANET) by setting certain combination of the actuators, we observe the effect on the
operation of the WDN. This effort sheds light on how the advanced algorithm in cyber
network can affect the purely water network through the integrated simulator and the
limitation of using EPANET to simulate the possible failures on the WDN. Furthermore, the
effort can validate the functionality that the game theory has in maintaining the equilibrium,
and how the equilibrium is reached in the EPANET reflected by the change of values in node
demand and flow level. The insight gained can be used to develop mitigation techniques that
harden the WDN against failures, ensuring a return on the considerable investment made in
adding cyberinfrastructure support to critical infrastructures.
Based on the completed work in the three parts, we conclude our contribution and present
our plan of research in the future.
2. Related work
As public safety concerns and prohibitive cost necessitate the use of modeling and simulation
for validation of intelligent environmental decision support systems (EDSSs), the utilization
of EDSSs in managing critical infrastructure has been investigated in numerous studies.
A general introduction to integrated decision support systems for environment planning
is provided in Kainuma et al. (1990). Applications of EDSSs include prevention of soil
salinization Xiao & Yimit (2008), regional environment risk management in municipal areas
Wang & Cheng (2010), and environmental degradation monitoring Simoes et al. (2003).
Examples particularly relevant to this book chapter are Xiao & Yimit (2008), which presents
an integrated EDSS for water resource utilization and groundwater control; and Serment
et al. (2006), which defines the major functionalities for an EDSS dedicated to the hydraulic
management of the Camargue ecosystem. Discussion on available models and tools, such
as GIS, and database management systems, is presented in Rennolls et al. (2004), which
also presents an application of biogeochemical modeling for sustainability management of
European forests.
Resource management algorithms have also been proposed for intelligent regulation. For
instance, hedging rules have been utilized to minimize the impact of drought by effectively
reducing the ongoing water supply to balance the target storage requirement Tu et al. (2003).
Applications of game theory include optimization of rate control in video coding Ahmad &
Luo (2006), allocation of power in frequency-selective unlicensed bands Xu et al. (2008), and
power control in communications MacKenzie & Wicker (2001). Most relevant to this book
108
4 Scientific and Engineering Applications Using MATLAB
Will-be-set-by-IN-TECH
chapter is the use of game theory in analyzing water resources for optimal allocation Yu-Peng
et al. (2006). Unlike our work, where the focus is to enable environmental management,
specifically water allocation, through the use of CPSs; the focus of Yu-Peng et al. (2006) is
on incorporating social and economic factors to provide a solution that maximizes the overall
value of water resources while satisfying both administrative resources allocation mandates
and consumer requirements.
This book chapter presents an EDSS, with the broader goal of applying the insights gained
to similar CPSs. Many CPSs, especially critical infrastructure systems, can be viewed as
commodity transport networks. WDNs are an example, as are smart grids and intelligent
transportation systems. The commodity transported varies from one domain to another, but
the systems share the goal of allocating limited resources under physical constraints, and
leverage the intelligent decision support provided by cyber infrastructure in achieving this
goal.
As an emerging research area, the body of literature specifically related to CPSs is limited. A
considerable fraction of related work examines critical infrastructure systems. The focus of
the majority of studies related to CPSs, e.g., Haimes & Jiang (2001); Pederson (2006); Rinaldi
(2004); Svendsen & Wolthusen (2007) is on interdependencies among different components
of critical infrastructure. A relatively comprehensive summary of modeling and simulation
techniques for critical infrastructure systems, an important category of CPSs, is provided
in Rinaldi (2004). Related challenges are enumerated in Pederson (2006), where system
complexity is identified as the main impediment to accurate characterization of CPSs. Other
challenges include the low probability of occurrence of critical events, differences in the time
scales associated with these events, and the difficulty of gathering data needed for accurate
modeling. Our work is one of few studies in the emerging field of CPSs to go beyond
qualitative characterization of the system to quantitative analysis.
Several challenges to the development of a generic framework for the design, modeling,
and simulation of CPSs are articulated in Kim & Mosse (2008). Features described as
desirable for such a framework include the integration of existing simulation tools, software
reusability, and graphical representation of the modeling and simulation environment. The
work presented in this book chapter meets all these criteria.
The study most closely related to the work presented in this book chapter is Al-Hammouri
et al. (2007), where a method is proposed for integration of the ns-2 network simulator with
the Modelica framework, a modeling language for large-scale physical systems. The paper
highlights the challenge of two-way synchronization of the simulators. The key difference
between this study and our work is that we link to a specialized simulator capable of
accurately representing the operation of the physical infrastructure, in this case a WDN, at
high resolution. The WDN simulator, and other related simulation tools are described in the
next section of this book chapter.
the end at higher hydraulic head to that at lower head, due to the effect of gravity. A negative
label for a flow indicates that its direction opposes that of the pipe.
In the WDN depicted in Fig. 2, the reservoir is providing water to the tank and a number
of different junctions. This topology can serve as a simple and abstract representation of a
lake that provides water to consuming entities spread throughout a city. The reservoir in this
figure always contributes water into the network, so its demand value is negative. The value
of the demand indicates the amount of water contributed, in this case 9884.69 gallons per
minute (GPM). The tank consumes the highest amount of water. Each junction is also labeled
with its demand value, and each pipe with its flow speed. The entire graph is color-coded to
simplify the categorization of demand or flow. The demand values of pumps and valves vary
in accordance with the nodes they control.
A more complex topology is depicted in Fig. 3, which shows a screen capture at hour 8:00 of a
24-hour simulation period. This figure also depicts node groupings, circled in green, that can
facilitate study of a subset of the nodes in the topology.
After simulating the system for the specified duration, EPANET can provide a report in graph,
table, or text form. Among the various reports available, the full report provides the most
comprehensive data, including the initial and updated values of all properties of the nodes
and links within each simulation time step (one hour by default). The water flow, pressure at
each node, depth of water in tanks and reservoirs, and concentration of chemical substances
can be tracked from the recorded data. Figs. 4 and 5 present snapshots of the link and node
information, respectively, of the full report.
simulation of the cyber layer of a WDN, as the decision support algorithms used are typically
implemented in a distributed fashion.
ns-2 USC Information Sciences Institute (2011), a public-domain discrete event simulator,
is the tentative choice for representing the communication network, an aspect of the cyber
infrastructure that is yet to be investigated.
3.3 Challenges in linking simulators for the cyber and physical networks
Accurate simulation of a CPS hinges on correctly recreating the information flow of Fig.1,
through the following iterative procedure:
112
8 Scientific and Engineering Applications Using MATLAB
Will-be-set-by-IN-TECH
0DWODE VLPXODWRUIRUF\EHULQIUDVWUXFWXUH
5XQGHFLVLRQ
3DUVHUHSRUW 2XWSXWWKHVH
VXSSRUWDOJRULWKPV
WRH[WUDFWLQSXW VHWWLQJVDVD
WRGHWHUPLQH
IRUDOJRULWKPV ,13ILOH
FRQWUROOHUVHWWLQJV
5XQ
6SHFLI\LQLWLDO 3URYLGHWKLV,13
(3$1(7DQG
:'1 ILOHWR(3$1(7DV
JHQHUDWHIXOO
FRQILJXUDWLRQ LQLWLDOFRQILJXUDWLRQ
UHSRUW
(3$1(7VLPXODWRUIRUSK\VLFDOLQIUDVWUXFWXUH
The first step in simulating an intelligent WDN is to specify the duration to be simulated
and the configuration of the physical infrastructure, e.g., topology and demand values, in
Integrated Cyber-Physical
Integrated Cyber-Physical Simulation
Simulation of Intelligent of Intelligent
Water Distribution Networks Water Distribution Networks 1139
EPANET. A 24-hour duration was selected for the simulation presented in this section. After
simulating the system for the specified duration, EPANET generates a full report that includes
information for all links and nodes for each time step (one hour by default), as shown in Figs.
4 and 5. The full report generated as the output file of EPANET is automatically saved as
a plain-text .NET file. This information includes values required as input by the decision
support algorithms of the cyber infrastructure, which in turn determine settings for physical
control elements such as valves.
To simulate the provision of sensor readings and other information about the physical
infrastructure to the cyber control system, the full report generated as output by EPANET
needs to be provided as input to Matlab. This necessitates pre-processing of the file, and
parsing of the data into the matrix form required by Matlab. A script using the textscan and
cell2mat commands can be defined within Matlab to carry out this pre-processing to generate
a separate matrix from the EPANET data for each entity (node or link) for each simulation
time step recorded in the full report, e.g., hour 1:00.
For simplicity, the simulation illustrated in this section was focused on node flow. The
controller (pump or valve) settings were determined by averaging the node demand within a
node group, which is a subset of nodes defined in EPANET. Fig. 3 shows a number of groups.
The same parsing approach can be used to extract additional data, e.g., water pressure or
concentration of a given chemical, from the EPANET report, as required by more sophisticated
decision support algorithms.
Each node group can reflect an associated group of consumers, such as residential nodes
in the south of a city. The only requirement is that each node group include at least one
controller (pump or valve), so controller settings determined by the cyber infrastructure can
be utilized in water allocation. The focus of the simulation in this section was integrated
simulation of the CPS, and as such, a simplistic approach was taken to water allocation, with
the goal of distributing the water as equitably as possible, subject to physical constraints on the
nodes. More intelligent decision support can be achieved through game-theoretic approaches
Yu-Peng Wang & Thian (2006), and it will be elaborated in Section 5.
Matlab generates a matrix of controller settings, which need to be provided to EPANET, as
they would be to the physical control elements in an actual WDN. A .INP file is required, in a
format identical to the original input provided to EPANET in the first step of the simulation,
with controller values updated to reflect the settings determined by the decision support
algorithm. A Matlab script utilizing the dlmwrite and fprintf commands can be used to generate
a .INP file with the format expected by EPANET.
In the final stage of the simulation, the .INP file generated by Matlab, which specifies settings
for various control elements, is used to initiate another execution of EPANET, closing the
physical-cyber-physical loop. The process can be repeated as necessary to simulate operation
of the WDN over multiple cycles of cyber control. Fig. 7 shows the file resulting from
execution of the water allocation algorithm for the node groups of Fig. 3. The result of
executing EPANET with the .INP file generated by Matlab is shown in Fig.8. As an example of
the manifestation of cyber control, the flow in the link connecting Junction1 (J1) and SOURCE,
marked with an arrow, has been reduced from 75-100 GPM (yellow) in Figure 3 to 50-75 GPM
(green) in Figure 8.
on the utilization of game theory for resource sharing and service provision in peer-to-peer
networks Gupta & Somani (2005).
If service is provided by player i in time period t, w is set to 1, otherwise 0. The reputation of all
players is initialized as 0 at time t = 0, and is defined as w at t = 1. Therefore, 0 ≤ R(t, i ) ≤ 1
is always maintained. In Equation 1, parameter a is a constant that captures the strength of
the “memory of the system,” i.e., the relative importance of current vs. past behavior of an
116
12 Scientific and Engineering Applications Using MATLAB
Will-be-set-by-IN-TECH
agent in determining its reputation. The notion of reputation is key in the game model, as it
affects the probability of receiving service for a player, and forms the incentive mechanism to
contribute service in the system. More detailed discussion is presented in Section 6.
The expected payoff value of electing to serve during time period t is defined as:
In Equation 2, the term (−C + R(t, Srv) ∗ U ) illustrates the tradeoff inherent to service
provision, namely, that cost of providing service as compared to the benefit of receiving
service. The term R(t, Srv) ∗ U reiterates that the probability of obtaining service in the current
time period depends on a player’s reputation. This payoff value of a player not only reflects
its current payoff after providing service, but also captures the potential to obtain service in
the next period, through the inclusion of R(t, Srv), which can be used as a health indicator
that reflects the capability of the player to gain service in the near future. When service is
provided, w = 1, and per Equation 1:
The equation reflects the “no contribution, no cost” case. When service is declined, w = 0,
and per Equation 1:
R(t, i ) = R(t − 1, i ) ∗ (1 − a) (5)
In a mixed-strategy Nash equilibrium of finite games, each player’s expected payoff should
be the same for all actions. In other words, the respective payoff values for {Srv} and {Dcln}
are equal:
Payoff(Srv) = Payoff( Dcln) (6)
Substituting from Equations 2 and 4 yields:
Incorporating the iterative definition of reputation, from Equations 3 and 5, the probability of
service provision, p, is determined as:
R ( t − 1) ∗ U (1 − a )
p= (8)
−C + 2R(t − 1) ∗ U (1 − a) + Ua
Several noteworthy points arise from the equations above. Firstly, p changes during each time
period, and is a function of the agent’s reputation at the end of the immediately preceding
period, R(t − 1). Secondly, recall that this is a mixed-strategy Nash equilibrium action profile,
where all players have the same p. Thirdly, we contend that this equilibrium is more stable
than the pure-strategy equilibrium discussed above, as self-interest will motivate agents to
eventually provide service in order to increase their chances of receiving service.
The agents are labeled Node i, Node j, and Node k, respectively. For each agent, the service
strategy is as shown in Table ??. The strategy shown in Table ?? does not exhaustively capture
all actions that could be taken by the three agents, but it provides a representative set of actions
over a non-trivial duration of ten time slots.
Time t Node i Node j Node k
1 Serve j Serve k Decline
2 Decline Serve i Decline
3 Serve k Decline Decline
4 Decline Decline Serve i
5 Serve k Decline Serve i
6 Serve j Decline Serve i
7 Serve j Serve i Decline
8 Decline Decline Decline
9 Decline Decline Serve j
10 Serve k Serve i Decline
According to the Table ??, we can summarize the strategy of each player, i, as Wi below:
• Wi = [1 0 1 0 1 1 1 0 0 1]
• Wj = [1 1 0 0 0 0 1 0 0 1]
• Wk = [0 0 0 1 1 1 0 0 1 0]
The configuration of initial values for the utility of obtaining service U and the cost of
providing service C is U/C = 80, with U = 800 and C = 10. The main reason to adopt
the ratio of utility to cost, U/C = 80, rather than their difference, U − C, is the normalization
inherent to use of the ratio. In civil engineering literature, water pricing has been approached
from a supply and demand perspective Brown & Rogers (2006); Cui-mei & Sui-qing (2009),
which is what U and C try to capture.
The U/C ratio can reflect whether the water resource is scarce or sufficient. U/C is low when
water is scarce, as serving a limited resource to other agents while maintaining sufficient
resources for own usage purpose will be expensive for an agent, leading to high C; and gaining
utility from other agents is difficult, leading to low U. Similarly, U/C is high when sufficient
water exists for all peer agents. Our initial choice of U/C = 80 for the simulation reflects a
non-draught situation. Simulation results for other values of U/C are presented in Lin et al.
(2011, to appear).
Integrated Cyber-Physical
Integrated Cyber-Physical Simulation
Simulation of Intelligent of Intelligent
Water Distribution Networks Water Distribution Networks 119
15
The main criteria for creating the topology include two aspects: the water distribution
network should have at least 3 actuators, either pump or valve, in charge of three different
areas, respectively; the water distribution network should have 3 reservoirs, representing
three agents to provide or retrieve water from their neighbors. Fig. 11 shows the grouped
nodes in the topology, which indicates what components are incorporated in the scope
managed by the particular agent. Each scope managed by one agent has one actuator.
to represent the state of serving water in one agent and 0 to represent the state of declining to
serve water (including retrieving water from other agents). Accordingly, in the first simulation
period, the script played by three agents is (1, 0, 0). Similarly as shown in the topology of
Fig. 9, we suppose the reservoir 1 is node i, and reservoir 8 and 9 are node j and node k,
respectively.
In terms of implementation, the water attributes (demand, pressure, head, flow, etc.) in
EPANET are controlled by the actuators (pump and valve). By sending the control command
to the actuator from the cyber infrastructure (implemented in MATLAB), we can configure the
serve/decline to serve operation of the node (reservoir). Because there are three actuators in
Fig. 11 with the open/close options, we can have totally 8 different combinations, shown in
Fig. 14.
The initial configuration (constraints) of the components (pump, valve, tank, node) can affect
the simulation result, and we set the initial values as following:
Integrated Cyber-Physical
Integrated Cyber-Physical Simulation
Simulation of Intelligent of Intelligent
Water Distribution Networks Water Distribution Networks 121
17
• All the three reservoirs (1, 8 and 9) have the total head of 1000 feet and elevation of 1000
feet.
• Valve 2 has 12 inch diameter, the type is PRV, loss coefficient is 0, and the fixed status is set
as none, as it can be open or closed.
• Pump 1 has pump curve 1 and the initial status is set as open.
• Valve 9 has 12 inch diameter, the type is PRV, loss coefficient is 0, and the fixed status is set
as none, as it can be open or closed.
• Tank 2 has elevation of 10 feet(the elevation above a common datum in feet of the bottom
shell of the tank) of 700 feet, initial level (the height of the water surface above the bottom
elevation of the tank at the start of the simulation), minimum level of 0 feet(the minimum
height in feet of the water surface above the bottom elevation that will be maintained;
the tank should not be allowed to drop below this level), maximum level of 20 feet (the
maximum height in feet of the water surface above the bottom elevation that will be
maintained; the tank should not be allowed to rise above this level) and 50 inch diameter.
• Junction 3 has elevation of 700 feet, base demand of 80 gpm, and its actual demand is
shown during simulation.
• Junction 4 has elevation of 500 feet, base demand of 75 gpm, and its actual demand is
shown during simulation.
• Junction 5 has elevation of 600 feet, base demand of 50 gpm, and its actual demand is
shown during simulation.
• Junction 6 has elevation of 500 feet, base demand of 20 gpm, and its actual demand is
shown during simulation.
• Junction 7 has elevation of 600 feet, base demand of 30 gpm, and its actual demand is
shown during simulation.
Subjected to the limited cases that actuators can manipulate the water flow and the constraints
of the capacity of pipe and node, such as the maximum flow the pipe can sustain for pump 1,
when we use the game theory to control the water resource on the EPANET, we need to take
these constraint factors into consideration and make decision accordingly.
Indicated by Fig. 14, we should avoid the failures generated by the two types of configuration
of the actuators. In another words, EPANET can not continue the simulation if pump 1, valve
2 and valve 9 are in the status of (open, open, open) or (open, close, open). This shows that the
command issued from the cyber simulator for controlling the actuators can lead to the errors
or malfunction of the underlying simulator for the physical network, and in this case, it is
EPANET.
According to Fig. 14, three patterns (strategy of the player) of water resource provision are
repeated consecutively, and they are (1, 0, 0), (1, 1, 1) and (1, 0, 1). We define the pattern as
serving pattern and the serving strategy similarly as Table ?? is the combination of the three
patterns. For example, if the initial serving pattern is (1, 0 ,0), the we configure the next serving
pattern as (1, 1, 1). There are multiple actuator setting methods to achieve this serving pattern,
for this case, we select the combination of (close, open, open), mapping with pump 1, valve
2 and valve9. All the rest of the configuration remains the same as initial configuration. The
generated control command file (input .INP file to EPANET) by MATLAB is shown in the
snapshot of as Fig. 15, which captures the part of actuator configuration. As shown in the
.INP file, the three actuators are configured as (close, open, open).
Integrated Cyber-Physical
Integrated Cyber-Physical Simulation
Simulation of Intelligent of Intelligent
Water Distribution Networks Water Distribution Networks 123
19
Fig. 17. Node demand (in GPM) at 0 hour when all three reservoirs are serving.
Fig. 18. Flow in the link (in GPM) at 0 hour when all three reservoirs are serving.
2. The reservoir in EPANET has the ability to provide infinite quantity of water, which could
be infeasible in real application case.
3. Although all the three reservoirs are providing water, the magnitude of provided water
quantity is different. Compared with reservoir 1, the water provided by reservoir 8 and
Integrated Cyber-Physical
Integrated Cyber-Physical Simulation
Simulation of Intelligent of Intelligent
Water Distribution Networks Water Distribution Networks 125
21
9 can almost be neglected, although at the beginning reservoir 8 and 9 are providing
more quantity of water. This change actually indicates the condition for reaching the
equilibrium in a water distribution network, i.e. the case that all the reservoirs (or players)
are consistently providing water is not an equilibrium or stable case, which conforms to
our previous analysis in subsection 5.2.1 on the pure equilibrium case.
4. The game theory in MATLAB is an supplemental intelligence onto the EPANET, and it is
an artificial manipulation for controlling the water rather than the hydraulic or physic
law. The purpose to define the reputation of player and the expected payoff value is
to investigate how the incentive mechanism for contributing service in the system can
affect the equilibrium in the water allocation. The more the player serves, the higher
reputation it can gain, and the higher probability it can gain water from other players.
The parameters in the game theory configure how the game will play among the players,
such as the probability that one player will serve in the next phase, but the strategy for
service game played among the players (i.e. which player serve and which player decline
to serve) determines the actual water allocation. In the combination of game theory and
the integrated CPS simulator, we directly use the strategy played among the players, and
set the configuration of actuators accordingly.
The effort of combining the CPS simulator and game theory shows the chain effect that the
advanced algorithm can issue a command of controlling the actuator, and the configuration of
the actuator can affect the simulation on the physical network. Sometimes the configuration
of the actuators may cause failures as indicated in Fig. 14, and this is due to the fact that
the physical components, such as pipes or tanks are subjected to the constraints which are
configured initially. The simulation reveals the risk that in real application, the calculated
configuration of the actuators can lead to the malfunctions of physical components, because
of the multiple constraints exerted on the components. This discovery can be used to
develop mitigation techniques that harden the WDN against failures, specifically, the design
of advanced computing algorithm on the cyber network needs to consider the multiple
constrains in the physical network, in order to ensure that adding the cyberinfrastructure to
support the operation of critical infrastructure will not bring serious reliability issues.
9. References
Ahmad, I. & Luo, J. (2006). On using game theory to optimize the rate control in video coding,
IEEE Transaction on Circuits and System for Video Technology 16.
Al-Hammouri, A., Liberatore, V., Al-Omari, H., Al-Qudah, Z., Branicky, M. S. & Agrawal,
D. (2007). A co-simulation platform for actuator networks, Proceedings of the 5th
International Conference on Embedded Networked Sensor Systems (SenSys ’07), ACM,
New York, NY, USA, pp. 383–384.
Brown, C. & Rogers, P. (2006). Effect of forecast-based pricing on irrigated agriculture: A
simulation, Journal of Water Resources Planning and Management 132(6): 403–413.
Cui-mei, L. & Sui-qing, L. (2009). Water price forecasting method based on marginal-cost
theory: a case study in China, World Environmental and Water Resources Congress 2009,
ASCE.
Dutch Ministry of Economics (2011). Waterspot.
URL: http://www.waterspot.nl/
Fudenberg, D. & Maskin, E. (1986). The Folk theorem in repeated games with discounting or
with incomplete information, Econometrica 54: 533–554.
Gupta, R. & Somani, A. K. (2005). Game theory as a tool to strategize as well as predict node’s
behavior in peer-to-peer networks, Proceedings of the 11th International Conference on
Parallel and Distributed Systems (ICPADS ’05).
Integrated Cyber-Physical
Integrated Cyber-Physical Simulation
Simulation of Intelligent of Intelligent
Water Distribution Networks Water Distribution Networks 127
23
United States Environmental Protection Agency (2011c). Ground Water and Rainmaker
Simulators.
URL: http://www.epa.state.oh.us/ddagw/SWEET/sweet simulators.html
United States Environmental Protection Agency (2011d). Water Quality Analysis Simulation
Program.
URL: http://www.epa.gov/athens/wwqtsc/html/wasp.html
USC Information Sciences Institute (2011). ns-2.
URL: http://nsnam.isi.edu/nsnam/index.php/User Information
Wang, B. & Cheng, H. (2010). Regional environmental risk management decision support
system based on optimization model for minhang district in shanghai, International
Conference on Challenges in Environmental Science and Computer Engineering.
Xiao, L. & Yimit, H. (2008). Environmental decision support system development for soil
salinization in the arid area oasis, Proceedings of the International Seminar on Business
and Information Management, IEEE Computer Society.
Xu, Y., Chen, W., Cao, Z. & Letaief, K. B. (2008). Game-theoretic analysis for power
allocation in frequency-selective unlicensed bands, Proceedings of the IEEE Global
Telecommunications Conference (GLOBECOM ’08).
Yu-Peng, W., Jian-Cang, X., Lintao, C., Ker, K. & Yew-Gan, T. (2006). Game analysis in
water resources optimal allocation, Proceedings of the International Conference on Hybrid
Information Technology (ICHIT ’06).
Yu-Peng Wang, Jian-Cang Xie, L. C. K. K. & Thian, Y.-G. (2006). Game analysis in water
resources optimal allocation, 2006 International Conference on Hybrid Information
Technology (ICHIT ’06).
8
1. Introduction
More recent technological advancements in microprocessor relays, combined with GPS
receivers for synchronization and accurate time stamping, is providing users advanced relay
systems with synchronized measurements, called synchrophasor measurements (IEEE
Power System Relaying Committee, 2002; Phadke, 2002; Marek, 2002). Synchrophasor
measurements together with advancements in digital communications, provides users with
the power system state at a rate of twenty times per second. Synchrophasor measurements
from different network locations when combined and processed in a central computer
system will provide users with the absolute phase angle difference between distant network
buses with an accuracy of tenths of an electrical degree. These types of central computer
systems, equipped with wide-area protection and control algorithms, will be able to better
address future system out-of-step conditions and other system problems because they will
have a better knowledge of what happens throughout the power system. In addition,
knowledge of online generation and load demand provided from synchrophasor
measurement systems will aid in balancing better the generation and load during islanding,
as well as minimizing load and generation shedding in order to preserve stability during
major system disturbances. Time synchronized phasor measurements provide a dynamic
view of a power system, combining these measurements in a central protection system
(CPS); this capability is used to set up a wide area control, protection and optimization
platform by means of new communication systems and (GPS), integrated application design
is shown in Figure 1. Figure 1 shows an integrated application design based on phasor
measuring units. When the system operates in extreme conditions, load shedding,
generation shedding, or system islanding must occur to prevent total system collapse
(Thorp et al., 1988; Centen et al., 1993; Guzman et al., 2002; Guzman et al., 2002). Typical
causes of system collapse are voltage instability or transient angle instability. These
instabilities can occur independently or jointly. In most cases, system wide-area disruptions
begin as a voltage stability problem. Because of a failure to take proper actions for the
system to recover, this voltage stability problem evolves into an angle stability problem.
New monitoring, protection, and communications technologies allow us to implement
economical local- and wide-area protection systems that minimize risk of wide-area system
disruptions or total system collapse.
130 Scientific and Engineering Applications Using MATLAB
control for better management of the system security through advanced control and
protection strategies.
The electricity supply industries need tools for dealing with system wide disturbances that
often cause cascading outages and widespread blackouts in power system networks. When
a major disturbance occurs, protection and control measures overtake the greatest role to
prevent further degradation of the system, restore the system to a normal state, and
minimize the impact of the disturbance. Electrical measurements of the system, which may
include synchronized phasors, are supplied to one or more wind farm controllers, which in
turn perform a control function improving the damping of electromechanical oscillations or
voltage performance in the utility system. The benefits are improved damping to
electromechanical oscillations and better voltage profile and ultimately more efficient
utilization of assets, reducing the necessity for installing new assets. Wide-area protection is
becoming an important issue and a challenging problem in the power industry (Wang et el.,
2005).
This study proposes a novel technique based on wide-area measurements for a power
system. The study is very vital and needed in the current state regarding the electrical utility
and the society as well to face future expansion of the electrical grid and to cover the
demand of the increasing growth and solving the problem of peak period. The study is very
beneficial also from the stability and security of the grid viewpoint in case of interconnection
with other countries.
This study presents a new approach for fault detection and classification for interconnected
system using the time synchronized phasor measurements. The scheme is depending on
comparing positive sequence voltage magnitudes for specified areas and positive sequence
current phase difference angles for each interconnected line between two areas on the
network. The chapter will cover all fault events for fault classification. The Matlab/simulink
program is extensively used to implement the idea. It uses to simulate the power system,
phase measurement unit function, synchronization process, fault detection and
classification.
centralization and decentralization. More and more dynamic functions are moving from
local and regional control centers toward central or national control centers. At the same
time we also observe more "intelligence" and "decision power" moving closer towards the
actual power system substations. Greater functional integration is being enclosed in
substation hardware. In view of global security of power systems, the action algorithms of
conventional backup protections possibly are not best choices because the operations of
individual relays are hardly coordinated each other. Therefore, the principle of the
protection design needs innovation to overcome the above problem (Yan et al., 2008;
Xiaoyang et al., 2008; Yangguang et al., 2010; Hui, et al., 2009).
Modem protection devices have sufficient computing and communications capabilities to
allow the implementation of many novel sophisticated protection principles. Therefore, a
novel wide-area backup protection system with fault classification is reported in this
chapter. This system is capable of acting as the substitution of conventional distributed
backup protections in substation. The architecture and algorithm of the system are also
introduced. To ensure the fast responsibility of such a system to the emergent events, the
communication requirements are discussed as well. Conclusively, the proposed system is
designed by two ways. First, in substation, concentrate some conventional backup
protection functions to an intelligent processing system; second, concentrate the coordinated
and optimized processing and controlling arithmetic of all backup protection in a region
into a regional processing unit. The communication of data among them is carried via optic-
fiber networks (Zhiyuan et al., 2009; 2009).
The proposed system comprises a master system and several local units. The system is
arranged as three layers. The bottom layer consists of PMUs with additional protection
functionality. The next layer consists of several Local Backup Protection Centers (LBPCs),
each of which interfaces directly with a number of PMUs. The top layer, System Backup
Protection Center (SBPC), acts as the coordinator for the LBPCs. Connected together via
fiber-optic communication links, these devices can process intelligent algorithms based on
data collected locally. The structure can be seen as Figure 2 (Seong et al., 2008).
Local part of the proposed wide-area backup system comprises PMUs and LBPC, which are
both installed in the station. PMUs are made up of DSP (Digital Signal Processor) and GPS
(Global Positioning System). The DSP measures instantaneous voltages and currents of
protected power system in real-time, and calculates the state variables, which provide vital
information for backup protection system. Then, power system variable data is transferred
to LBPC, LBPC samples digital inputs, and pre-processes the analog and digital signals, and
then deliver pre processed results to SBPC via fiber-optic communications between the
substation and SBPC. SBPC will be installed at Regional Control Center, and integrates
various well-developed functions, such as data acquisition via communication, system
monitoring, fault location, security analyzing, making tripping strategy and descending
strategy to LBPCs, also, SBPC can do post-event playback and post-event data analyses.
These functions employ PMUs to fulfil real-time demand (Testa et al., 2004).
LBPCs perform the correlative operations on the spot when receiving the strategy from
SBPC, then fault will be isolated rationally. While it is not possible to prevent all
contingencies that may lead to power system collapse, a wide area backup protection system
that provides a reliable security isolated scheme and optimized coordinated actions is able
to mitigate or prevent large area disturbances.
A Novel Wide Area Protection Classification
Technique for Interconnected Power Grids Based on MATLAB Simulation 133
Background SBPC
SCADA
Center Processor
Optical Fiber
Network
Local
part
LBPC LBPC LBPC
3.2 Inputs to control and protection systems (Moxley et al., 2007; Phadke et al., 2009;
Jetti et al., 2006)
The state of the power system is represented by several network parameters. Thresholds,
trends, patterns, and sudden changes of these parameters provide key information to detect
an emergency. Some of the key system parameters which constitute the possible inputs to
improved protection and control systems are:
Active power flows in the network - If the limits on active power are violated, the system is
in a viability crisis. For the overloaded transformer, a loss-of-life occurs. Thus guidance for
loading is established to assure a long life. The limit for the transmission line loading is set
by transient and steady-state stability conditions (usually long lines), voltage collapse
conditions (usually medium lines), and thermal conditions (usually short lines).
Voltage magnitude and reactive power flows - The voltages in the power system as well as
sudden voltage changes need to be contained within a small range. The voltage and
reactive power and their rate-of-change can provide valuable information on voltage
instability.
Angles between buses - Stability limits for every line will be satisfied, if the difference in
angles across the line does not exceed a certain limit. Detection of the out-of-step condition
can prevent instability, and, consequently, cascading.
Impedance - Unstable swing, stable swing, and fault condition may be detected and
distinguished by observing behavior of the impedance loci at the local bus. A typical out-of-
step blocking or tripping scheme is accomplished by "blinders" or circles in R-X diagram
and timers.
Resistance and rate-of-change of resistance - These parameters may be used to speed-up the
out-of-step detection.
Frequency - Frequency deviation from the nominal value is a result of power imbalance. In
modern interconnected systems, frequency deviation usually occurs in the islanded area (a
definite indicator of "in extremis" crisis).
Rate of change of frequency - Unlike frequency, rate of change of frequency is an
instantaneous indicator of power deficiency in the islanded area. The oscillatory nature of
the rate of change of frequency needs to be considered in utilizing this feature.
A Novel Wide Area Protection Classification
Technique for Interconnected Power Grids Based on MATLAB Simulation 135
Spinning reserve - The spinning reserve quantity, distribution, and the speed of its'
dynamical response are factors that influence the effectiveness of the spinning reserve
during an emergency. The speed of the dynamic response for the hydro units the first few
seconds after a demand is made is relatively slow compared to thermal units.
Consequently, the spinning reserve needs to be distributed throughout the system on both
hydro and thermal units. The spinning reserve needs to be considered in load shedding
schemes to optimize shed load.
Load - Load is a non-linear function of voltage and frequency. These changes in load impact
power system imbalance and frequency behavior. Further, load changes with the season and
the time of the day. In addition, underfrequency load shedding programs specify percent of
the total load that should be shed at each step. As load changes, actual load for shedding
does not correspond to planned load.
Relays and breaker status - Operation of the protective relays (desired or undesired) and
network configuration have an essential impact to disturbance propagation. If undesired
operation may be avoided by detecting hidden failures or by adapting relay settings to
prevailing system conditions, unwanted transition of the system to a less desirable
emergency state may be prevented. Further, equipment unavailability because of
maintenance and testing needs to be recognized and considered.
Modeling of the power network is required to simulate disturbances and to choose features
that will be extracted. The disturbance in the power network usually develops gradually;
however some phenomena, such as a rise of transient instability, can develop in a fraction of
second. Selection of appropriate power network analysis tools is important (load flow,
transient stability, mid and long term dynamic models, EMTP, etc.).
modern information technology can contribute to improve this situation. Our proposal to
review the present protection strategy to counteract large area disturbances addresses the
potentials that are derived from the advances in system operational, protection and control
techniques. It will be explained how the application of numerical technology can avoid
catastrophic disturbances to occur or at least to keep the impact of single fault within certain
limits.
In contrast to the requirements for protection relays designed to protect individual plant
objects, system protection schemes intended to prevent voltage or frequency instabilities
have to cope with the loss of generation on a large scale and/or loss of one or more
transmission lines. Information technology offers digital applications, in terms of numerical
adaptive protection relays, integrated disturbance recorders and fast broadband
communication much greater functionality and overall efficiencies than conventional
analogue techniques.
TL
CT
VT
transform, the signal is assumed to exist from time - to + but in the DFT, the signal exists for
a small duration of time (called window). The components of different frequencies
determined by the DFT analysis can be combined to recreate the original waveform.
incorporated into any power system design or analysis, as excess delays could ruin any
control procedures adopted to stabilize the power grid.
Although more and more control systems are being implemented in a distributed fashion
with networked communication, the unavoidable time delays in such systems impact the
achievable performance. Delays due to the use of PMUs and the communication link
involved are due primarily to the following reasons:
Transducer delays: Voltage transducers (VT) and current transducers (CT) are used to
measure the RMS voltages and currents respectively, at the instant of sampling.
Window size of the DFT: Window size of the DFT is the number of samples required to
compute the phasors using DFT.
Processing time: The processing time is the time required in converting the transducer data
into phasor information with the help of DFT.
Data size of the PMU output: Data size of the PMU message is the size of the information
bits contained in the data frame, header frame and the configuration frame.
Multiplexing and transitions: Transitions between the communication link and the data
processing equipment leads to delays that are caused at the instances when data is retrieved
or emitted by the communication link.
Communication link involved: The type of communication link and the physical distance
involved in transmitting the PMU output to the central processing unit can add to the delay.
Data concentrators: Data concentrators are primarily data collecting centers located at the
central processing unit and are responsible for collecting all the PMU data that is
transmitted over the communication link.
9. Conventional problems
The distance relays which are widely applied in the protection today and involve the
determination of impedance achieve operating times of the order of a period of the power
system frequency. A distance relay is designed to only operate for faults occurring between
the relay location and the selected reach point, and remains stable for all faults outside this
region or zone (Horowitz et al., 2009).
The resistance of the fault arc takes the fault impedance outside the relay’s tripping
characteristic and, hence, it does not detect this condition. Alternatively, it is only picked up
either by zone 2 or zone 3 in which case tripping will be unacceptably delayed (Eissa, 2009).
The distance relays are based on standalone decision, while each relay operates
independently according to three different zone of operation, see Figure 7.
The mal-operation or fail-to trip of protection is determined as one of the origins to raise and
propagate major power system disturbances (Tang, et al., 2006). A vast majority of relay
mal-operations is unwanted trips and have been shown to propagate major disturbances.
Backup protections in fault clearance system have the task to operate only when the primary
protection fails to operate or when the primary protection is temporarily out of service. The
recent complexity and enlargement of power systems makes it difficult to coordinate
operation times and reaches among relays. In the areas of power system automation and
substation automation, there are two different trends: centralization and decentralization.
More and more dynamic functions are moving from local and regional control centers
toward central or national control centers. At the same time we also observe more
“intelligence” and “decision power” moving closer towards the actual power system
substations. Greater functional integration is being enclosed in substation hardware. In view
of global security of power systems, the action algorithms of conventional backup
protections possibly are not best choices because the operations of individual relays are
hardly coordinated each other. Therefore, the principle of the protection design needs
innovation to overcome the above problem. Modern protection devices have sufficient
computing and communications capabilities to allow the implementation of many novel
sophisticated protection principles. Therefore, a novel wide-area backup protection system
is reported in this paper.
This system is capable of acting as the substitution of conventional distributed backup
protections in substation. To ensure the fast responsibility of such a system to the emergent
events, the communication requirements are discussed as well. Conclusively, the proposed
system is designed by two ways. First, in substation, concentrate some conventional backup
protection functions to an intelligent processing system; second, concentrate the coordinated
and optimized processing and controlling arithmetic of all backup protection in a region
into a regional processing unit.
The communication of data among them is carried via optic-fiber networks. The relay
decision is based on collected and shared data through communication network. The
suggested technique satisfies high degree of reliability and stability while it is based on
shared decision rather than stand alone decision. The suggested technique can see all the
power system area and can deal with the transmission lines as unit protection, see Figure 8.
The primary purpose of these systems is to improve disturbance monitoring and system
event analysis. These measurements have been sited to monitor large generating sites, major
transmission paths, and significant control points. Synchronized Phasor measurements
provide all significant state measurements including voltage magnitude, voltage phase
angle, and frequency.
bus for each area. This can result in the minimum voltage value that indicates the nearest
area to the fault. In addition to that, the absolute differences of the positive sequence current
angles are calculated for all lines connected with the faulted area. These absolute angles are
compared to each other. The maximum absolute angle difference value is selected to
identify the faulted line. The above two keys of operation can be mathematically described
as follows:
ф12
ф21
ф13
ф31
ф23
ф32
ф35
ф53
ф24
ф42
ф45
ф54
Fig. 11. The Matlab/Simulink block diagram responsible of the selection of the faulted line.
A Novel Wide Area Protection Classification
Technique for Interconnected Power Grids Based on MATLAB Simulation 145
Part of the 500/220 kV Egyptian interconnected electrical network is used for the study; five
main buses that represent five different areas with 500 kV are selected to verify the
suggested technique. Figure 14 shows the selected five areas from the overall network. In
the single line diagram, each bus represents the selected area in the simulation that can
connect the 500 kV network with 220 kV network through three single phase 500/220 kV
power transformers. The system is simulated using the Matlab/Simulink with a sampling
frequency of 20 kHz for a system operating at a frequency of 50 Hz. By means of measuring
positive sequence magnitude of three phase to ground voltage and positive sequence
current angle difference between sending and receiving ends, we make a new criteria to deal
with faults (single, double and three phase to ground faults). This new criteria will detect
faults and select the nearest area to the fault. Also, the new criteria will distinguish between
internal and external faults in the interconnected lines.
Fig. 12. Maximum absolute difference of positive sequence current angle due to internal
fault.
Fig. 13. Maximum absolute difference of positive sequance current angle due to external
fault.
146 Scientific and Engineering Applications Using MATLAB
Shoubra
Area 2 Korimat
Motamadia
Maghagh
a
A.zaabal
Demo
Samalu
t
Area 5
Selected
W.Houf
Asuit
C.South
Area 4
Area 1 500 kV
Tebbin South
220 kV
Fig. 14. The selected five areas from the overall network.
voltage signals measured from Kurimat bus bar is recorded and displayed in Figure 16. The
3-phase current signals for all transmission lines away from Kurimat are recorded and
displayed in Figure 17.
1
λ zero = I a + Ib + Ic (5)
3
The fault is classified as a ground fault if λzero is greater than threshold value. Once the fault
has been classified, the specific fault type is determined by comparing all rate of change of
phase currents with an predetermined disturbance detection pickup, The fault typing
148 Scientific and Engineering Applications Using MATLAB
software module is employed when the measured parameter exceeds a prescribed threshold
(e.g., when a measured correct exceeds an overcurrent pickup).
Fig. 17. Three phase current signals for all lines connected to Kurimat.
Fig. 19. Positive sequence current angle absolute differences for all lines connected to the
faulted area.
Discrete Fourier transform block computes the fundamental value of the input phase
current signal over a running window of one cycle of the specified fundamental frequency
as shown in Figure 20. First and second outputs return respectively the magnitude and
phase degrees of the fundamental. The magnitude is taken as a percentage from its steady
state value. The rate of change of this percentage is compared with a threshold value. For
the first cycle of simulation, the outputs are held constant to the value specified by the
parameter "Initial input".
As shown in Figure 21, the input three-phase current is used to calculate zero sequence
components to classify the fault type. Then each phase current signal is taken as a
percentage from its steady state value. Then the rate of change of the percentage of phase
current magnitude is compared by a threshold value to identify the faulted phase.
Fig. 22. Fault recorder at each terminal display that fault type is three phase fault.
13. Conclusion
The chapter outlines a novel idea for fault detection and classification using Phasor
measurement units in a wide area system. The idea has successfully identified the faulted
line on a large power interconnected system. The idea descried in this paper represents a
new state-of-art in the field of interconnected grid protection and classification. The idea is
based on sharing data from many PMUs. The new idea also calssified the fault types for the
interconnected system. The idea used a center protection unit for collecting the data and
issued the tripping signal. The idea is implemented and investigated using the powerful
Matlab/Simulink package. Power system configuration, fault detection, fault calculation,
discrimination, and classification are achieved through the Matlab/Simulink program.
14. Acknowledgment
This chapter is written, revised and analyzed by the author. Some parts of the chapter are
cited from MSc dissertation done in Helwan University-Faculty of Engineering-Cairo-Egypt
by Eng. Mohamed Magdy. The MSc dissertation is supervised by Prof. Mohammed El
Shahat Masoud and the author of the chapter.
15. References
Phadke, A. (2002). Synchronized Phasor Measurements a historical over view. Virginia
Polytechnic Institute and American Electric Power Bhargava.0-7803-7525-4/02/IEEE
Marek Zima. (2002). Special Protection Schemes in Electric Power Systems. Literature survey,
(June 2002).
Thorp, J. S; Phadke A. G; Horowitz, S. W. & Begovic, M. M. (1988). Some Applications of
Phasor Measurements to Adaptive Protection. IEEE Transactions on Power Systems,
Vol. 3, No. 2, (May 1988).
Centeno, V; De La Ree, J; Phadke, A. G.; Mitchell, G; Murphy, J; Burnett, R; (1993). Adaptive
Out-of-Step Relaying Using Phasor Measurement Techniques. IEEE Computer
Applications in Power, (October 1993).
Guzman, A & Tziouvaras, D (2002). Local and Wide-Area Network Protection Systems
improve power system reliabili, Schweitzer Engineering Laboratories, E. O. ty. Inc.
Pullman, WA, USA, Ken E. Martin Bonneville Power Administration Vancouver, WA,
USA, Aviable from, www.naspi.org/resources/archive/prtt/waps_wprc04.pdf
Guzman, A; Mooney, J. B; Benmouyal, G; Fischer, N. (2002). Transmission Line Protection
System for Increasing Power System Requirements. 55th Annual Conference for
Protective Relay Engineers, College Station, Texas, (April 2002).
Gjerde, P. O & Mangelrod, R. (2001). Optimisation of Protection Performance During System
Disturbances. CIGRE 2001 SC 34 Colloquium, (Septemper 2001).
Larsson, M & Rehtanz, C. (2002). Predictive Frequency Stability Control Based on Wide-
Area Phasor Measurements. Proc. 2002 IEEE Power Engineering Society Summer
Meeting, (2002).
Zima, M; Korba, M & Andersson, G. (2003). Power Systems Voltage Emergency Control
Approach Using Trajectory Sensitivities. 2003 IEEE Conference on Control
Application, CCA, Istanbul, (June 2003).
152 Scientific and Engineering Applications Using MATLAB
Rehtanz C & Westermann, D. (2002). Wide Area Measurement and Control System for
Increasing Transmission Capacity in Deregulated Energy Markets. Proc. 2002 Power
Systems Computation Conference in Sevilla, (2002).
Rehtanz, C. (2001). Wide Area Protection System and Online Stability Assessment based on
Phasor Measurements. Bulk Power System Dynamics and Control - V, Onomichi
City, Japan, (Augest 2001).
Zhang, N & Kezunovic, M. Verifying the Protection System Operation Using An Advanced
Fault Analysis Tool Combined with the Event Tree Analysis. Northern American
Power Symposium, NAPS 2004 Moscow, Idaho, (August 2004).
Lin, Y. L; Liu, C. W. & Chen, C. S. (2004). A new PMU-based fault detection/location
technique for transmission lines with consideration of arcing fault discrimination—
Part I: Theory and algorithms. IEEE Trans. Power Del., Vol. 19, No. 4, (October
2004), pp. 1587–1593.
Yu, C. S; Liu, C. W; Yu, S. L. & Jiang, J. A. (2002). A new PMU-based fault location algorithm
for series compensated lines. IEEE Transaction on Power Deivery, Jan. 2002, Vol. 17,
No. 1, pp. 33–46.
Wang, et el. (2005). Design of a novel wide-area backup protection system, Proceedings of
IEEE/PES Transmission Distribution Conference Asia Pac. Dalian, China, (2005), pp. 1-6.
Yan Wang; Yan-Xia Zhang; Song-Xiao Xu; (2008). Wide-area protection against chain over-
load trip based on multi-agent technology, Interntional Conference on Machine
Learning and Cybernetics, pp. 1548 – 1552, Volume: 3, Digital Object Identifier:
10.1109/ICMLC.2008.4620652, 2008.
Xiaoyang Tong; Xiaoru Wang; Li Ding (2008). Study of information model for wide-area
backup protection agent in substation based on IEC61850, Third International
Conference on Electric Utility Deregulation and Restructuring and Power Technologies,
pp. 2212 – 2216, Digital Object Identifier: 10.1109/DRPT.2008.4523778
Yangguang Wang; Xianggen Yin; Dahai You; (2010). Agent-based wide area protection with
high fault tolerance, International Conference on the Modelling, Identification and
Control (ICMIC), (2010), pp. 739 – 744
Hui Sun; Qianjin Liu (2009). Research of the Structure and Working Mechanisms of Wide-
Area Backup Protection Agent. 1st International Conference on Information Science
and Engineering (ICISE), (2009), Digital Object Identifier: 10.1109/ICISE.2009.926,
(2009), pp. 5037 – 5042.
Zhiyuan Duan; Zhang, C.; Hu, Z.; Sun, Y. (2009). Design of Robust Controller of
Interconnected Power System Considering Signals Transmission Delay.
International Workshop on Intelligent Systems and Applications, (2009). ISA 2009.
Digital Object Identifier: 10.1109/IWISA.2009.5072823, (2009), pp. 1 - 5
Zhiyuan Duan; Chengxue Zhang; Zhijian Hu; Yuanyuan Zhang (2009). Robust Control of
Interconnected Power System Based on WAMS Considering Signals Transmission
Delay, Power and Energy Engineering Conference, (2009). APPEEC 2009, Asia-Pacific,
Digital Object Identifier: 10.1109/APPEEC.2009.4918729, (2009), pp. 1 – 4
Seong-Jeong Rim; Seong-Il Lim; Seung-Jae Lee (2008). Multi-agent based reliability
enhancement scheme for IEC61850 substation,
Transmission and Distribution Conference and Exposition, (2008), T&D. IEEE/PES,
Digital Object Identifier: 10.1109/TDC.2008.4517294, (2008), pp. 1 – 6
A Novel Wide Area Protection Classification
Technique for Interconnected Power Grids Based on MATLAB Simulation 153
Testa, S. & Chou, W. (2004). The distributed data center: front-end solutions, IT Professional
Conference, Vol. 6, No. 3, Digital Object Identifier: 10.1109/MITP.2004.24, (2004),
pp. 26 – 32
Begovic, M. (2004). Wide area protection and emergency control. Power Systems Conference
and Exposition, (2004). IEEE PES Digital Object Identifier:
10.1109/PSCE.2004.1397488, (2004), pp. 1776 - 1777
Terzija, V.; Cai, D.; Valverde, G.; Regulski, P.; Vaccaro, A.; Osborne, M. & Fitch, J. (2010).
Flexible Wide Area Monitoring, Protection and Control applications in future
power networks Developments in Power System Protection (DPSP 2010), 10th IET
International Conference on Managing the Change, Digital Object Identifier:
10.1049/cp.2010.0361 (2010), pp. 1 – 5
Terzija, V.; Valverde, G.; Deyu Cai; Regulski, P.; Madani, V.; Fitch, J.; Skok, S.; Begovic, M.
M. & Phadke, A. (2011). Wide-Area Monitoring, Protection, and Control of Future
Electric Power Networks, Proceedings of the IEEE, Vol. 99, No. 1, Digital Object
Identifier: 10.1109/JPROC.2010.2060450, (2011), pp. 80 - 93
Moxley, R. & Wronski, M. (2007). Using time error differential measurement in protection
applications, Power Systems Conference: Advanced Metering, Protection, Control,
Communication, and Distributed Resources, 2007. PSC 2007 Digital Object Identifier:
10.1109/PSAMP.2007.4740918, (2007), pp. 278 - 283
Phadke, A. G. & Kasztenny, B. (2009). Synchronized Phasor and Frequency Measurement
Under Transient Conditions, IEEE Transactions on Power Delivery, Vol. 24, No. 1,
Digital Object Identifier: 10.1109/TPWRD.2008.2002665, (2009), pp. 89 - 95
Jetti, S. R. & Venayagamoorthy, G. K. (2006). Real-Time Implementation of a Dual Function
Neuron based Wide Area SVC Damping Controller, Industry Applications
Conference, (2006), 41st IAS Annual Meeting. Conference Record of the 2006 IEEE,
Vol. 2, Digital Object Identifier: 10.1109/IAS.2006.256598, (2006), pp. 672 - 678
Yi Hu & Novosel, D. (2006). Challenges in Implementing a Large-Scale PMU System Power
System Technology, International Conference on Power Con 2006. Digital Object
Identifier: 10.1109/ICPST.2006.321829, (2006), pp. 1 - 7
Wenxin Liu; Cartes, D. A. & Venayagamoorthy, G. K. (2006). Particle Swarm Optimization
based Defensive Islanding of Large Scale Power System, International Joint
Conference on Neural Networks, IJCNN, (2006), Digital Object Identifier:
10.1109/IJCNN.2006.246642
(2006), pp. 1719 - 1725
Xiupeng Guan; Yuanzhang Sun & Lin Cheng (2007). A load parameter identification method
based on wide area measurement, International Power Engineering Conference, 2007.
IPEC 2007, (2007), pp. 431 - 436
Zhang, G.; Lee, S.; Carroll, R.; Jian Zuo; Beard, L. & Yilu Liu (2010). Wide area power system
visualization using real-time synchrophasor measurements, Power and Energy
Society General Meeting, 2010 IEEE Digital Object Identifier:
10.1109/PES.2010.5588188, (2010), pp. 1 - 7
Moraes, R.; Volskis, H. & Yi Hu (2008). Deploying a large-scale PMU system for the
Brazilian interconnected power system, Third International Conference on Electric
Utility Deregulation and Restructuring and Power Technologies, (2008), DRPT
2008, Digital Object Identifier: 10.1109/DRPT.2008.4523392, (2008), pp. 143 - 149
154 Scientific and Engineering Applications Using MATLAB
Sybille, G. & Hoang Le-Huy (2000). Digital simulation of power systems and power
electronics using the MATLAB/Simulink Power System Blockset, IEEE Power
Engineering Society Winter Meeting, (2000), Vol. 4 Digital Object Identifier:
10.1109/PESW.2000.847358, (2000), pp. 2973- 2981
De La Ree, J.; Centeno, V.; Thorp, J.S. & Phadke, A. G. (2010). Synchronized Phasor
Measurement Applications in Power Systems, IEEE Transactions on Smart Grid, Vol.
1, No. 1, Digital Object Identifier: 10.1109/TSG.2010.2044815, (2010), pp. 20 - 27
Hall, I.; Beaumont, P.G.; Baber, G.P.; Shuto, I.; Saga, M.; Okuno, K. & Uo, H. (2003), New line
current differential relay using GPS synchronization, Power Tech Conference
Proceedings, 2003 IEEE Bologna, Vol. 3 Digital Object Identifier:
10.1109/PTC.2003.1304464, 2003
Naduvathuparambil, B; Valenti, M. C. & Feliachi, A. (2002). Communication delays in wide
area measurement systems, Lane Dept. of Comp. Sci. & Elect. Eng., West Virginia
University, Morgantown, WV, 2002. Avialble from
www.csee.wvu.edu/~mvalenti/documents/SSST_02.pdf
Klump R. & Wilson, R. E. (2005). Visualizing real-time security threats using hybrid
SCADA/PMU measurement displays, Proc. 38th Hawaii Int. Conf. System Sci., 2005,
p. 55c.
Bhargava, B. & Salazar, A. (2008). Use of Synchronized Phasor Measurement system for
monitoring power system stability and system dynamics in real-time, Power and
Energy Society General Meeting-Conversion and Delivery of Electrical Energy in
the 21st Century, 2008 IEEE Digital Object Identifier: 10.1109/PES.2008.4596963,
2008, pp. 1 - 8
Horowitz S. H. & Phake A. G. (2009). Power System Relaying, Taunton Somerset, U.K.
Research Studies Press
Eissa, M. M. (2009). New principle for transmission line protection using phase portrait
plane. IET Generation Transmission Distribution, (2009), Vol. 3, No. 1, pp. 49–56
Tang J. & McLaren, P. G. (2006). A wide area differential backup protection scheme for
Shipboard application, IEEE Transaction on Power Delivery, (2006), Vol. 21, No. 3,
pp. 1121-1129
Eissa, M. M. (2008). Development and investigation of a new high-speed directional relay
using field data, IEEE Transaction on Power Delivery, (2008), Vol. 23, No.3, pp. 1302–
1309.
Eissa, M. M. (2009). A new digital feed circuit protection using directional element, IEEE
Transaction on Power Delivery, (2009), Vol. 24, No. 2, pp. 531–537
Eissa, M. M. (2005). Evaluation of a new current directional protection technique using field
data, IEEE Transaction on Power Delivery, (2005), Vol. 20, No. 2, pp. 566–572
Protection of Interconnected Electrical Networks Using Phasor Synchronized Measuring
Technique, Ph.D. dissertation, Helwan University-Faculty of Engineering at
Helwan, Cairo, Egypt, (2008)
9
1. Introduction
The need of ultrawideband (UWB) antennas with omni-directional coverage is increasing
for both military and commercial applications (Wiesbeck et. al., 2009; Minin, 2010). The UWB
radio technology promises high resolution radar applications, sensor networks with a large
number of sensors for industrial or home surveillance as well as high data-rate
communication over short range for personal area networks. With a need for antennas with
the characteristics of broad bandwidth and small electrical size, conical antenna structures
have been a focus of research because of its broad bandwidth and omni-directional radiation
pattern (Maloney & Smith, 1993; Sandler & King, 1994; Yu & Li, 2008; Palud et. al., 2008). The
bi-conical antenna exhibits a very stable omni-directional radiation pattern in the plane
normal to the dipole axis together with an excellent transient response. However, the
feeding with a usual coaxial cable requires a balun, which transforms the asymmetric mode
of the feed line into a symmetric mode at the feed point. For the coaxial balun the ultra wide
bandwidth demands very high precision in the manufacturing process in order to get a
good and stable matching especially for the high frequencies. The mono-cone antenna as
asymmetric structure does not need any balun for an asymmetric feed line but it needs an
infinite ground plane, which in reality can only be approximated. The theory of wide-angle
conical antennas has been developed sufficiently to permit calculation of the transfer
functions relating source voltage to radiated field and incident field to load voltage over the
range of frequencies required in the study of transients. Such calculations were
demonstrated in (Harrison & Williams, 1965).
Due to their three-dimensional configurations, conical antennas are bulky and difficult to
fabricate, integrate, and reconfigure. Moreover, since conventional conical antennas
comprise of free-standing metal, they are typically heavy in order to achieve sufficient
mechanical stability. Several configurations have been proposed to improve conical
antennas’ mechanical performance (Ma et. al., 2009; Zhou et. al., 2009, Kliros et. al., 2010a).
Resistive loading for conical antennas, which is investigated in (Maloney 1993), does not
constitute the optimal solution as it reduces the antennas’ efficiency. Recently, investigations
have been carried out on configurations that employ a dielectric or magnetic material to
cover the conical antenna (Gentili et. al., 2004, Lu, 2007). Dielectric and magnetic coating of
156 Scientific and Engineering Applications Using MATLAB
the radiating cone, enables making the antenna electrically smaller and more rugged while
maintaining a wide band input impedance.
In this chapter, we present a Finite Difference Time Domain (FDTD) code in spherical
coordinates implemented in MATLAB in order to simulate the performance of dielectric
covered conical antennas. MATLAB provides an interactive environment for algorithm
development, data post-processing and visualization. The spherical FDTD equations can be
found using a modified Yee cell in spherical coordinates (Fusco, 1990). Spherical Berenger’s
perfectrly matched layer (PML) is applied as absorbing boundary condition where a
parabolic conductivity profile in the spherical PML-region is used (Berenger, 1996). A
unique feature of the PML is that electromagnetic waves of arbitrary incidence, polarization
and frequency are matched at the boundary in a reflectionless manner. Results concerning
time evolution of the radiated electromagnetic field, the return loss, input impedance,
maximum gain as well as far-field radiation patterns across an extended bandwidth, are
presented. A time domain study has also been performed to characterize the antenna’s
behaviour in case an UWB pulse is used. For evaluating waveform distortions caused by the
antenna, we examine the degree of similarity between source pulse and received pulse
waveforms in several propagation directions. The effect of the dielectric spherical cover on
the antenna’s performance is investigated. The author mostly worked in MATLAB version
7.4 and the related sample codes are provided in the Appendix.
∂Ε
∇×Η = ε ⋅ +σ ⋅Ε (2)
∂t
where ε and μ are the permittivity and permeability respectively and, σ and σ* the electric
and magnetic conductivity of the propagation space respectively. These two vector
equations are the general equations governing the antenna operation. Because the conical
antenna is a revolutionary symmetric structure, the three dimensional problem can be
reduced to two-dimensional problem. In spherical coordinates, due to rotational symmetry,
Eq. (1) and (2) lead to the following three scalar partial differential equations:
∂Ηφ σ* ∂ ∂
( sin θ ⋅ Εr )
1
=− ⋅ Ηφ − sin θ ( r ⋅ Εθ ) − (3)
∂t μ μ ⋅ r ⋅ sin θ ∂r ∂θ
∂Εr σ 1 ∂
∂t
= − ⋅ Εr +
ε ε ⋅ r ⋅ sin θ ∂θ
sin θ ⋅ Ηφ ( ) (4)
∂Εθ σ 1 ∂
∂t
= − ⋅ Εθ −
ε ε ⋅ r ∂r
r ⋅ Ηφ ((5) )
To obtain a discrete set of the continuous differential equations, the central difference
approximation is used on both the time and space first-order partial derivatives. The entire
computational space is a collection of modified Yee unit cells (Yee, 1966). In our modified
Yee’s scheme, the computational space is subdivided by using an orthogonal mesh in
spherical coordinates. The electric fields are located along the edges of the cells, while the
magnetic fields are positioned at the centers of these cells. Using the well-known half time
step notation in all locations and after some rearrangements, a set of finite difference field
forms for Eqs. (3)-(5) follows (Kliros et. al., 2010a, 2010b):
Cb (i ) sin ( ( j + 1 / 2) Δθ ) n + 1 n+
1
Ern + 1 ( i , j ) = C a (i ) Ern ( i , j ) + ⋅ Hφ 2 ( i , j ) − Hφ 2 ( i , j − 1 ) (6)
( i + 1 / 2 ) Δθ sin (( j − 1 / 2) Δθ )
n+ 1 1
i + 1 / 2 n+ 2
Εθn + 1 (i , j ) = C a (i ) Εθn (i , j ) + C b (i ) ⋅ Hφ 2 (i − 1, j ) − φ ( i , j )
⋅ H (7)
i − 1/2
n+
1
n−
1
Db (i ) sin ( ( j + 1 ) Δθ ) n
Ηφ 2
( i , j ) = Da (i ) Ηφ 2 ( i , j ) + ⋅ Εr ( i , j + 1 ) − Εrn ( i , j )
( i + 1 / 2 ) Δθ sin ( j Δθ )
(8)
i + 1 n
Eθ ( i + 1, j ) − Eθ ( i , j )
n
− Db ( i )
i
where
σ ( i ) Δt Δt
1−
2ε (i ) ε (i ) Δr
C a (i ) = , Cb (i ) = (9)
σ (i ) Δt σ (i ) Δt
1+ 1+
2ε (i ) 2ε (i )
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 159
σ * (i )Δt Δt
1−
2 μ(i ) μ (i ) Δr
Da (i ) = , Db (i ) = (10)
σ * (i )Δt σ ∗ (i )Δt
1+ 1+
2 μ(i) 2 μ (i )
and Δr , Δθ represent the step size in the r- and θ- directions, respectively. Superscript n
signifies that the quantities are to be evaluated at t = nΔt, and, i and j represent the point
(iΔr, jΔθ) in the spherical grid. The half time steps indicate that the fields E and H are
calculated alternately. The maximum time step is limited by the stability Courant’s criterion
(Fusco, 1990):
r ( Δr )( Δθ )
Δt ≤ min (11)
c ( Δr )2 + ( r Δθ )2
where c is the velocity of the light in free space.
Fig. 2. Spherical perfectly matched layer at the edge of the simulation space.
160 Scientific and Engineering Applications Using MATLAB
∂Ηφθ 1 ∗ 1 ∂ ( Εr sin θ )
= −σ θ Ηφθ + (12)
∂t μ r sin θ ∂θ
∂Ηφ r 1 ∗ 1 ∂ ( r ⋅ Εθ )
=− σ r Ηφ r + (13)
∂t μ r ∂r
∂Εr 1
= −σ θ Ε r +
1 ∂ Ηφ r + Ηφθ sin θ (( ) ) (14)
∂t ε r sin θ ∂θ
∂Εθ 1
= − σ r Εθ +
1 ∂ r Ηφ r + Ηφθ (( )) (15)
∂t ε r ∂r
ii. Create spherical FDTD equations from the above revised Maxwell equations:
n+
1
n−
1
Dbθ (i ) sin ( ( j + 1 ) Δθ ) n
Ηφθ 2 ( i , j ) = Daθ (i ) Ηφθ 2 ( i , j ) + ⋅ Εr ( i , j + 1 ) − Εnr ( i , j ) (16)
( i + 1 / 2 ) Δθ sin ( )
j Δ θ
1 1
i+1 n
( i , j ) = Dar ( i ) Η φ r 2 ( i , j ) − Dbr ( i )
n+ n−
Eθ ( i + 1, j ) − Eθ ( i , j )
n
Ηφ r 2 (17)
i
C bθ ( i ) sin ( ( j + 1 / 2) Δ θ ) n + 1 n+
1
Ern + 1 ( i , j ) = Cαθ ( i ) Ern ( i , j ) + ⋅ H φ 2 ( i , j ) − H φ 2 ( i , j − 1 ) (18)
( i + 1 / 2 ) Δ θ sin (( j − 1 / 2) Δθ )
n+ 1 i + 1/2 n+
1
Εθn + 1 ( i , j ) = C ar ( i ) Εθn ( i , j ) + C br ( i ) ⋅ H φ 2 ( i − 1, j ) − ⋅ H φ ( i , j )
2 (19)
i − 1 / 2
where
σ r Δt Δt
1− ε Δr
Cα r = 2 ε , C br = (20)
σ Δt σ Δt
1+ r 1+ r
2ε 2ε
σ r∗ Δt Δt
1−
2μ μ Δr (21)
Dar = , Dbr =
σ ∗ Δt σ ∗ Δt
1+ r 1+ r
2μ 2μ
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 161
σ θ Δt Δt
1−
C aθ = 2 ε , C = εΔ θ (22)
bθ
σ Δt σ Δt
1+ θ 1+ θ
2ε 2ε
σ θ∗ Δt Δt
1−
2μ μΔ θ (23)
Daθ = , Dbθ =
σ θ∗ Δt σ ∗ Δt
1+ 1+ θ
2μ 2μ
and Ηφ = Ηφ r + Ηφθ in the last two equations.
iii. for a given number N of PMLs, calculation of free-space conductivities, σ 0 and σ 0* , the
final conductivities σ N and σ N* and the conductivity profile of each PML.
According to (Berenger, 1996), for a desired conductivity profile σ (r ) of thickness δ, the
reflection factor at normal incidence R(0) is given by:
2 δ
R(0) = exp − σ (r ) dr (24)
ε c 0
and consequently, the reflection factor for a wave at arbitrary incidence, θ, is
ε 0c ln ( R ( 0 ) )
σ0 = − (26)
2 4 Δr N PML
3
2
i
(27)
σ ( i ) = σ max , σ max = 24σ 0 N PML
δ
The correct conductivity profile is calculated automatically in our code for any desired
reflection factor and number of perfectly matched layers.
Vs (t ) − I in (t ) Rs
Eθn (i , j ) = − (28)
b sin(θ )ln(b/a)
162 Scientific and Engineering Applications Using MATLAB
where
n − 1/2
Vs (n Δt ) − Rs I in (i , j )
Eθn (i , j ) = − (30)
b ⋅ sin( jΔθ )ln(b/a)
n − 1/2
I in (i , j ) = 2π (i Δr )sin ( ( j + 1 / 2 ) Δθ ) ⋅ Hφn − 1 2 ( i + 1 / 2, j + 1 / 2 ) (31)
The above field is a spherical source extending from the inner conductor of the coaxial line
to the outer conductor. Consequently, as the voltage Vs(t) steps forward in time, it drives the
base of the conical antenna with the above spherical field. Then, the antenna radiates the
resulting wave following the time evolution described by Eq. (6)-(8).
Vin ( f )
exp( − jπ f Δt )
Zin ( f ) = (32)
I in ( f )
where the exponential term accounts for the half-time step difference between the electric
and magnetic field computation.
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 163
The results of input impedance are then used to obtain the return loss characteristics of the
antenna. Thus, the Return Loss S11 (dB) of the antenna is given by
Zin ( f ) − Zline
Γ( f ) = (34)
Zin ( f ) + Zline
is the frequency dependent reflection coefficient. From the calculated reflection coefficient,
the VSWR can be calculated as follows:
1 + Γ( f )
VSWR = (35)
1 − Γ( f )
The bandwidth of the antenna is the frequency range corresponding to a reflection
coefficient of the antenna less than or equal to 1/3 that leads to VSWR ≤ 2.
To calculate the conical antenna gain, the far electric field in the desired direction must be
determined as a function of frequency. Since the electric far-field is computed so that the 1/r
amplitude factor and the propagation delay are suppressed, the antenna gain relative to a
lossless isotropic antenna in θm – direction is given by
2
1 E( f ,θ m )
G( f ,θ m ) = (36)
2η Pin / 4π
where E( f ,θ m ) is the peak value of the Fourier transform of the pulsed far field radiated in
the θm-direction, η the characteristic space impedance and Pin the steady-state input power
at each frequency given by
1
Pin ( f ) = ∗
Re Vin ( f )I in ( f ) (37)
2
Directivity curves can be also computed. Directivity of an antenna is defined as the ratio of
the radiation intensity in a given direction from the antenna to the radiation intensity
averaged over all directions. The directivity in θm-direction is given by the ratio of the
integral of Poynting vector with the value of electric field E( f ,θ m ) to the actual value of the
integral:
2
2 E( f , θ m )
D( f ,θ m ) = (38)
π
E ( f ,θ )E( f ,θ )sin θ dθ
∗
The FDTD cell dimensions are Δr=3 mm and Δθ=1ο. The antenna sits on top of a perfectly
conducting ground plane that extends 360o in all directions for a distance of Rm=10 . Just
before the maximum radial distance Rm is reached, the simulation space is terminated by a
PML section of thickness 20Δr. The maximum reflection coefficient at normal incidence is
chosen to be R(0) = 10-14. The time step is taken Δt=0.2 psec, sufficient to satisfy Courant’s
criterion. An UWB Gaussian pulse (FWHM = 64 psec) modulated by a continuous sine wave
carrier of frequency fc is used in our simulations, that is,
t − t0 )
2
Vs (t ) = Vm exp − ( sin ( 2π f c ( t − t0 ) ) (39)
td2
where τ = 64 p sec, t0 = 4 × τ , f c = 6.5 GHz and Vm = 0.1V . The UWB excitation pulse driving
the conical antenna is depicted in Fig.5.
1
0.8
0.6
0.4
no rm aliz e d V s (t)
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
t (nsec)
Fig. 5. Excitation UWB pulse driving the conical antenna.
the frequency bandwidth. The received pulses are a bit larger due to the fact that the
antenna has filtered all frequencies outside the impedance bandwidth. The longer duration
of the received pulses indicates lower achievable data-rates while the shape distortion can
make the detection process more difficult. Proper channel models can be used to study these
effects (Molisch, 2003).
(a)
(b)
(c)
Fig. 6. FDTD simulation snapshots of the electric field strength after 2000 time steps (a), 5000
time steps (b) and 9000 time steps (c).
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 167
Fig. 7. Excitation and radiated pulses versus the elevation angle. It can clearly be seen that
the radiated signals are elevation angle dependent.
200
o
180 θ 0 =25
o
θ 0 =35
160
o
θ 0 =47
Re(Z in ) ( Ω )
140 o
θ 0 =55
120
100
80
60
40
2 4 6 8 10 12 14 16 18
frequency (GHz)
Figure 10, shows the evolution of maximum gain in the elevation plane (E-plane) versus
frequency. Gain gradually increases with frequency from 10 dBi to about 14 dBi in the
frequency range from 5.5 to 8.5 GHz and remains almost frequency independent at 14 dBi in
the frequency range from 8.5 to 17 GHz.
Furthermore, the power radiation patterns in the elevation plane are calculated in the above
frequency range, although for brevity, only the patterns at 4.5, 6.5, 8 and 10 GHz are shown
in Fig.11. Obviously, the power radiation patterns present quasi-perfect omni-directional
(monopole-like) behaviour but gradually degrade with increasing frequency. The radiation
lobe enlarges downwards up to 6.5 GHz, and above 10 GHz a ‘null’ appears near θ=35o
while the beamwidth decreases slowly with frequency. These variations are attributed to the
fact that the antenna’s electrical size increase with frequency. It is also observed that the
radiation patterns are slightly upward looking. This feature could be useful for radar sensor
network applications.
4
o
θ=25
o
3.5 θ=35
o
θ=47
o
θ=55
3
VSWR (dB)
2.5
1.5
1
2 4 6 8 10 12 14 16 18
frequency (GHz)
20
15
Maximum Gain (dBi)
10
0
o
θ =25
-5 o
θ =35
o
-10 θ =47
o
θ =55
-15
-20
2 4 6 8 10 12 14 16 18
frequency (GHz)
Fig. 10. Maximum gain for various flare angles and εr=3.0.
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 169
120
ε r =1.0
110
ε r =2.2
100 ε r =3.0
ε r =4.4
Impedance ( Ω )
90
ε r =9.8
80
70
60
50
40
30
2 4 6 8 10 12 14 16 18
frequency (GHz)
Fig. 12. Input impedance for various flare angles and εr=3.0.
As it is seen in Figs.12 and 13, the ultra wide-band characteristics of the antenna are not
sensitive to the variation of dielectric constant εr. Some ripples appeared in both real part of
170 Scientific and Engineering Applications Using MATLAB
impedance and VSWR are smoothed out when the dielectric cover is present. As it is seen,
when εr is reasonably small, the antenna’s input impedance remains close to a constant (~50
Ω) within a wide frequency band. However, because of the dielectric-air interface, the
reflection and scattering at the end of the conical radiator is stronger making the antenna
less matched to free space. Consequently, a wide range of dielectric materials can be used to
construct the spherical cover of the antenna. Obviously, low conductivity material is
preferred in order to minimize the dielectric loss. Figure 14 illustrates the frequency
dependence of the maximum gain for covers of different dielectric constant. As it is seen, the
maximum gain remains almost frequency independent in a wide range of frequency for all
materials. However, as the dielectric constant increases, the maximum gain decreases.
2.5
ε r =1.0
ε r =2.2
ε r =3.0
ε r =4.4
2
ε r =9.8
VSWR
1.5
1
2 4 6 8 10 12 14 16 18
frequency (GHz)
Fig. 13. Simulated VSWR for various flare angles and εr=3.0.
20
15
10
Maximum Gain (dBi)
ε r =1.0
-5
ε r =2.2
-10 ε r =3.0
ε r =4.4
-15 ε r =9.8
-20
2 4 6 8 10 12 14 16 18
frequency (GHz)
Fig. 14. Maximum gain for various flare angles and εr=3.0.
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 171
4. Conclusion
In this Chapter, we present a FDTD code in spherical coordinates implemented in MATLAB
in order to simulate the radiation characteristics of conical antennas. MATLAB provides an
interactive environment for algorithm development, data post-processing and visualization.
The spherical FDTD equations can be found using a modified Yee cell in spherical
coordinates. Spherical Berenger’s perfectrly matched layer (PML) is applied as absorbing
boundary condition where a parabolic conductivity profile in the spherical PML-region is
used.
The code is used to design and simulate a conical antenna covered by a spherical dielectric
structure and placed above a large ground plane. This quasi-planar antenna is mechanically
stable and, relative easy to build and integrate with the planar circuits. Parametric studies
lead to the optimum values of cone’s arm length =45 mm and flare angle θ0=470 for 50 Ω
matched impedance. This design achieves an impedance bandwidth from 5.5 to 17 GHz,
with stable radiation patterns over this bandwidth. The radiation patterns are monopole-like
and their frequency dependence is small in the whole UWB frequency band. A time domain
study has shown that the antenna distorts the excitation pulse in a moderate way.
It is observed that, the ultra wide-band characteristics of the antenna are not sensitive to the
variation of dielectric constant εr of the spherical cover. Consequently, a wide range of
dielectric materials can be used to construct the spherical cover of the antenna.
Our study suggests that a spherical dielectric covered conical antenna holds sufficient
potential as a low-profile antenna with very wideband characteristics. A very important
need is to verify more of the simulation results with experimental measurements which will
be reported in a future communication.
5. Acknowledgment
The author wishes to thank his former graduate students G. Kyritsis and D. Touzloudis
(now First-Lieutenants at Hellenic Air-Force) for improving and testing the first version of
FDTD-code as well as some post-processing MATLAB-codes during their diploma thesis.
6. Appendix
===================================================================
tmax=16384; % we use 2^{14} steps!
k=0; mu_0 = 4*pi*1e-7; eps_0 = 1e-9/(36*pi);
impedance_free=120*pi; c = 3e+8; Zline=50;
% ==================== Define the antenna dimensions =======================
ant_length= 15; % ant_length in units (dr) (mm).
ant_angle = 48; % flare angle:'ant_angle'Add (+1)for Matlab
e_rel=3.0; % relative dielectric constant of substrate
tan_loss=0.0009d0; % tangent loss of dielectric substrate
sub_length=ant_length; % defines dielectric substrate length
e_sub=e_rel*eps_0;
impedance_medium=impedance_free/sqrt(e_rel);
radial_view=121; % give the position (r,theta) to view the far-fields
theta_view=46;
% =================== Definitions and Constants ============================
dr = 0.003; % radial step
172 Scientific and Engineering Applications Using MATLAB
else
sigma(I,J) = sigmaPML(I + Npml - radial_max);
sigmaM(I,J) = sigmaMPML(I + Npml - radial_max);
end
end
if ((I <=sub_length) & (J > ant_angle )& (J< (theta_max-1)))
sigma(I,J)= sigma_epsilon;
sigmaM(I,J)= sigmaM_epsilon;
end
end
end
for I = 1:radial_max %initialize, zero the fields for all free nodes
for J= 1:theta_max
if ((I >= ant_length)|(J >= ant_angle))
Er(I,J) = 0; Et(I,J) = 0; Hp(I,J) = 0;
end
end
end
g1 = dt/(2*mu_0); g2 = dt/(2*eps_0);
g3 = dt/(dr*eps_0); g4 = dt/(dr*mu_0);
I_factor = 2*pi*b;
V_factor = (log(sin(ant_angle*pi/180))-log(1-cos(ant_angle*pi/180)))/log(2);
%============================ Begin of time iterations =====================
t = 0;
fid=fopen('Etime.dat','wt');
while (t < (tmax*dt))
t = t + 0.5*dt; % For the first half time-step, update the H-field.
for I = 3:(radial_max-Npml)
for J = 2:(theta_max-1)
if ((I >= ant_length)|(J >= ant_angle))
g5 = sin(J*dth)/sin((J-1)*dth);
Da = (1-sigmaM(I,J)*g1)/(1+sigmaM(I,J)*g1);
Db = g4/(1+sigmaM(I,J)*g1);
ER1 = (g5*Er(I,J+1)-Er(I,J))/((I-1/2)*dth);
ET1 = (I/(I-1))*Et(I+1,J) - Et(I,J);
Hp(I,J) = Da*Hp(I,J) + Db*(ER1 - ET1);
end
end
end
J = theta_max; % H-field at the ground plane
for I = 3:(radial_max-Npml)
Da = (1-sigmaM(I,J)*g1)/(1+sigmaM(I,J)*g1);
Db = g4/(1+sigmaM(I,J)*g1);
ET1 = (I/(I-1))*Et(I+1,J) - Et(I,J);
Hp(I,J) = Da*Hp(I,J)- Db*ET1;
end
for I = ant_length:(radial_max-Npml)
Hp(I,1) = Hp(I,2); % H-field along the line of symmetry
end
layer = 1; % Update the H-field in the PML region
for I = (radial_max-Npml+1):radial_max
Dar = (1-sigmaMPML(layer)*g1)/(1+sigmaMPML(layer)*g1);
Dbr = g4/(1+sigmaMPML(layer)*g1);
Dat = Dar; Dbt = Dbr;
for J = 2:(theta_max-1)
g5 = sin(J*dth)/sin((J-1)*dth);
ER1 = (g5*Er(I,J+1)-Er(I,J))/((I-1/2)*dth);
174 Scientific and Engineering Applications Using MATLAB
if(I == radial_max)
ET1 = 0;
else
ET1 = (I/(I-1))*Et(I+1,J)- Et(I,J);
end
Hpr(layer,J) = Dar*Hpr(layer,J) - Dbr*ET1;
Hpt(layer,J) = Dat*Hpt(layer,J) + Dbt*ER1;
Hp(I,J) = Hpr(layer,J) + Hpt(layer,J);
end
layer = layer + 1;
end
J = theta_max; % Update H-field in the PML along the ground plane
layer = 1;
for I = (radial_max-Npml+1):radial_max
Dar = (1-sigmaMPML(layer)*g1)/(1+sigmaMPML(layer)*g1);
Dbr = g4/(1+sigmaMPML(layer)*g1);
Dat = Dar;
if(I == radial_max)
ET1 = 0;
else
ET1 = (I/(I-1))*Et(I+1,J) - Et(I,J);
end
Hpr(layer,J) = Dar*Hpr(layer,J) - Dbr*ET1;
Hpt(layer,J) = Dat*Hpt(layer,J);
Hp(I,J) = Hpr(layer,J) + Hpt(layer,J);
layer = layer + 1;
end
t = t + 0.5*dt; % For the second half time-step, update the E-fields.
Vsource= Vmax*exp(-((t-t0)/td)^2)*sin(2*pi*fc*(t-t0));
Vin = 0; Iins = 0;
for J = ant_angle:theta_max
Iin = I_factor*sin((J-1/2)*dth)*Hp(3,J);
Vdrv = Vsource- Rs*Iin;
if ((t >= tsample) & (J > ant_angle))
Vin = Vin + b*Et(3,J)*dth;
Iins = Iins + Iin;
if(J == (theta_max-1))
Iins = Iins/(theta_max - ant_angle -1);
exVin = Vin; exIin = Iin;
end
end
Et(3,J) = Vdrv*(1/(b*log(2)))/sin((J-1/2)*dth);
end
k=k+1;
piVin(k)=Vin;
piIin(k)=Iin;
for I = 4:(radial_max-Npml) % Step E-fields for free nodes
for J = 2:(theta_max-1)
if ((I >= ant_length)|(J >= ant_angle))
g6 = sin((J-1/2)*dth)/sin((J-3/2)*dth);
Ca = (1 - g2*sigma(I,J))/(1 + g2*sigma(I,J));
Cb = g3/(1 + g2*sigma(I,J));
HPHI1 = (g6*Hp(I,J)- Hp(I,J-1))/((I-1/2)*dth);
HPHI2 = Hp(I-1,J)-((I-1/2)/(I-3/2))*Hp(I,J);
Er(I,J) = Ca*Er(I,J) + Cb*HPHI1;
Et(I,J) = Ca*Et(I,J) + Cb*HPHI2;
end
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 175
end
end
% ================ Introduction of spherical dielectric cover ==============
for I = 4:sub_length
for J = ant_angle:theta_max-1
g6 = sin((J-0.5)*dth)/sin((J-1.5)*dth);
g2_sub = 1/(1*e_sub/eps_0);
g3_sub = 1/(1*e_sub/eps_0);
Ca = (1 - g2_sub*sigma(I,J))/(1+ g2_sub*sigma(I,J));
Cb = g3_sub/(1 + g2_sub*sigma(I,J));
HPHI1 = (g6*Hp(I,J)- Hp(I,J-1))/((I-0.5)*dth);
HPHI2 = Hp(I-1,J)-((I-0.5)/(I-1.5))*Hp(I,J);
Er(I,J) = Ca*Er(I,J) + Cb*HPHI1;
Et(I,J) = Ca*Et(I,J) + Cb*HPHI2;
end
end
J = theta_max;E % Compute the fields at the ground plane
for I = 4:(radial_max-Npml)
Ca = (1 - g2*sigma(I,J))/(1+ g2*sigma(I,J));
Cb = g3/(1 + g2*sigma(I,J));
HPHI2 = Hp(I-1,J) - ((I-1/2)/(I-3/2))*Hp(I,J);
Er(I,J) = 0;
Et(I,J) = Ca*Et(I,J) + Cb*HPHI2;
end
J = theta_max; % Compute the E-fields for the substrate domain
for I = 4:sub_length
g2_sub = 1/(1*e_sub/eps_0);
g3_sub = 1/(1*e_sub/eps_0);
Ca = (1 - g2_sub*sigma(I,J))/(1+ g2_sub*sigma(I,J));
Cb = g3_sub/(1 + g2_sub*sigma(I,J));
HPHI2 = Hp(I-1,J) - ((I-1/2)/(I-3/2))*Hp(I,J);
Er(I,J) = 0;
Et(I,J) = Ca*Et(I,J) + Cb*HPHI2;
end
for I = ant_length:(radial_max-Npml) %E-field along line of symmetry
Ca = (1 - g2*sigma(I,1))/(1 + g2*sigma(I,1));
Cb = g3/(1+g2*sigma(I,1));
HPHI1 = (Hp(I,2) - Hp(I,1))/((I-1/2)*dth);
Er(I,1) = Ca*Er(I,1) + Cb*HPHI1;
Et(I,1) = 0;
end
layer = 1; % Update the E-fields in the PML region
for I = (radial_max-Npml+1):radial_max
for J = 2:(theta_max-1)
g6 = sin((J-1/2)*dth)/sin((J-3/2)*dth);
Car = (1 - g2*sigmaPML(layer))/(1 + g2*sigmaPML(layer));
Cbr = g3/(1 + g2*sigmaPML(layer));
Cat = Car; Cbt = Cbr;
HPHI1 = (g6*Hp(I,J) - Hp(I,J-1))/((I-1/2)*dth);
HPHI2 = (Hp(I-1,J) - ((I-1/2)/(I-3/2))*Hp(I,J));
Er(I,J) = Cat*Er(I,J) + Cbt*HPHI1;
Et(I,J) = Car*Et(I,J) + Cbr*HPHI2;
end
layer = layer + 1;
end
layer = 1; %E-field in the PML region at the ground plane
J = theta_max;
176 Scientific and Engineering Applications Using MATLAB
for I = (radial_max-Npml+1):radial_max
Car = (1 - g2*sigmaPML(layer))/(1 + g2*sigmaPML(layer));
Cbr = g3/(1 + g2*sigmaPML(layer));
HPHI2 = (Hp(I-1,J) - ((I-1/2)/(I-3/2))*Hp(I,J));
Er(I,J) = 0; Et(I,J) = Car*Et(I,J) + Cbr*HPHI2;
layer = layer + 1;
end
Etheta(k)=Et(radial_view,theta_view);
Radial_view=151; tt=t/1e-12;
for sJ=1:theta_max-1
EthetaT(sJ,k)=Et(radial_view,sJ);
end
end
fclose(fid);
table=[tt; Et(radial_view, theta_view)];
%========== Spherical to rectangular coordinates transformation ============
for I=1:radial_max
for J=1:theta_max-1
x= radial_max +round((I*sin((J-1)*dth)));
x2 = (2*radial_max+1) - x; y = 1 + round((I*cos((J-1)*dth)));
Ecart(x,y) = Et(I,J); Ecart(x2,y) = Et(I,J);
end
end
for I=1:(2*radial_max-1)
for J = 1:radial_max
EcartNew(I,J) = Ecart(I,J);
end
end
Imin = 2; Imax = 2*radial_max - 2; Jmin = 2; Jmax = radial_max-1;
for I = Imin:Imax
for J = Jmin:Jmax
if ((Ecart(I,J)=0)&(radial_max*cos((I-radial_max)*dth*91/151)+25>= J))
ItempLo = I-1; ItempHi = I+1;
while ((Ecart(ItempLo,J) == 0) & (ItempLo > 1))
ItempLo = ItempLo-1;
end
while ((Ecart(ItempHi,J) == 0) & (ItempHi < 2*radial_max -1))
ItempHi = ItempHi + 1;
end
M = Ecart(ItempLo,J); N = Ecart(ItempHi,J);
if(M == 0)
temp1 = N;
elseif(N == 0)
temp1 = M;
else
temp1 = sign(M+N)*sqrt(abs(M*N));
end
JtempLo = J-1; JtempHi = J+1;
while ((Ecart(I,JtempLo)= 0) & (JtempLo > 1))
JtempLo = JtempLo-1;
end
while ((Ecart(I,JtempHi)= 0) & (JtempHi < radial_max))
JtempHi = JtempHi + 1;
end
M = Ecart(I,JtempLo); N = Ecart(I,JtempHi);
if(M == 0)
temp2 = N;
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 177
elseif(N == 0)
temp2 = M;
else
temp2 = sign(M+N)*sqrt(abs(M*N));
end
if (temp1==0)
EcartNew(I,J) = temp2;
elseif(temp2 == 0)
EcartNew(I,J) = temp1;
else
EcartNew(I,J) = sign(temp1+temp2)*sqrt(abs(temp1*temp2));
end
end
end
EcartNew(1,radial_max)=1.0; EcartNew(2,radial_max)= -1.0;
end
%======================== 3D time evolution of E-field =====================
I=1:2*radial_max-1; J=1:radial_max;
x(I) = I; y(J) = J;
surfl(y(J), x(I), EcartNew(I,J))
xlabel('Y-axis cm'), ylabel('X-axis cm')
zlabel('Etheta V/m')
shading interp; colormap bone
%========================== Frequency-Domain Analysis =====================
t_val=0.2e-12:0.2e-12:0.2e-12*tmax;
dt = t_val(2) - t_val(1); t0 = t_val(1); V_val=piVin; I_val=piIin;
iplot = 1; ipad = 0; omega_plot_max =2*pi*10e+9; omin=2; omax=61;
central_f=14; % omega (14)corresponds to frequency 6.5 GHz!
Vft_power,Vft,omega,iflag] =
get_Fourier_transform(dt,t0,V_val,iplot,ipad,omega_plot_max);%Call FFT
iflag=0; N=length(Vft);
[Ift_power,Ift,omega,iflag] = ...
get_Fourier_transform(dt,t0,I_val,iplot,ipad,omega_plot_max);
Z_FT=complex(zeros(1,N),zeros(1,N)); cutoff=1e-25;
for m=1:N
if(abs(Ift(m)< cutoff))
Z_FT(m)=0;
else
Z_FT(m)=Vft(m)/Ift(m);
end
end
j=1:N; Z=Z_FT(j); % Impedance Matrix
imath=sqrt(-1); phase=exp(-imath*0.5*omega*dt);
ZTR1=Vft./Ift; ZTR=ZTR1.*phase;
Re=real(ZTR); Im=imag(ZTR); Zabs=abs(ZTR);
plot(omega(omin:omax)/(2*pi*1e9),Re(omin:omax),'k-.')
xlabel('frequency (GHz)'); ylabel('Impedance (Ohms)');
GAMMA=(ZTR-Zline)./(ZTR+Zline); % Reflection coefficient
RLoss=20*log10(abs(GAMMA)); % Input Return Loss in (dB)
VSWR=(1+abs(GAMMA))./(1-abs(GAMMA)); % Voltage standing wave ratio
plot(omega(omin:omax)/(2*pi*1e9),VSWR(omin:omax),'k-.')
xlabel('frequency (GHz)'); ylabel('VSWR');
FeedPower=0.5*real(Vft.*conj(Ift)); % Input Power versus frequency
omeg=omega(central_f);
plot(omega(omin:omax)/(2*pi*1e9),RLoss(omin:omax),'k-.')
xlabel('Frequency (GHz)'); ylabel('Return Loss (dB)');
178 Scientific and Engineering Applications Using MATLAB
7. References
Berenger J.P. (1996). Perfectly Matched Layer for the FDTD Solution of Wave-Structure
Interaction Problems. IEEE Transactions on Antennas &Propagation, Vol.44, No. 1, pp.
110-117
Brocato, R. W. (2004). FDTD simulation tools for UWB antenna analysis, Sandia Report,
SAND2004-6577, Sandia National Laboratories.
Fusco M. (1990). FDTD Algorithm in Curvilinear Coordinates. IEEE Transactions on Antennas
& Propagation, 38, pp. 76-89
Harrison C. W., Jr. & C. S. Williams, Jr. (1965). Transients in Wide-Angle Conical Antennas,
IEEE Transactions on. Antennas & Propagation, Vol. AP-13, pp. 236-246
Simulated Performance of Conical Antennas
Using Matlab-Based Finite-Difference Time Domain (FDTD) Code 179
Gentili G. B., Cerretelli M. and Cecchi L. (2004). Coated conical antennas for automotive
application. Journal Electromagnetic Waves & Applications, Vol.18, No.1, pp. 85-97
Katz, D.S., Thiele, E.T & Taflove, A. (1994). Validation and extension to three dimensions of
the Berenger PML absorbing boundary condition for FDTD meshes. IEEE
Microwave & Guided Wave Letters, Vol, 4, No.8, pp. 268 - 270
Kliros G. S., Kyritsis G. & Touzloudis D. (2010), Dielectric-EBG covered conical antenna for
UWB applications. COMPEL: Internationlal Journal for Computation & Mathematics in
Electrical and Electronic Engineering, Vol. 29, no. 4, pp. 1134-1143
Kliros G. S., Kyritsis G. & Touzloudis D. (2010). FDTD Analysis of a Conical Microstrip
Antenna on EBG–substrate. Journal of Applied Electromagnetism (JAE), Vol. 12, No.1,
pp.1-8, ISSN: 1109-1606
Liang X. & Wah M. C. Y. (2000). Low-profile broadband omni-directional monopole
antenna. Microwave & Optical Technology Letters, Vol. 25, No. 2, pp.135–138
Liu G. & Grimes C. (1999). Spherical-coordinate FDTD analysis of conical antennas mounted
above finite ground planes. Microwave & Optical Technology Letters, Vol.23, pp.78-82
Lu M., Bredow J., Jung S. & Tjuatja S. (2007). A Quasi-Planar Wide Band Conical Antenna.
In: Ultra-Wide Band Short-Pulse Electromagnetics, Baum, C. E., Stone, A. P., Tyo, J. S.
(Eds.), pp. 25-32, Springer, ISBN 978-0-387-73045-5, New York
Ma J., Yin Y.Z., Zhou S.G. & Zhao L.Y. (2009). Design of a new wideband low-profile conical
antenna. Microwave & Optical Technology Letters, Vol.51, No.11, pp. 2620–2623
Maloney J. G. & Smith G. S. (1993). Optimization of a conical antenna for pulse radiation: an
efficient design using resistive loading. IEEE Transactions on Antennas &
Propagation, Vol. 41, No.7, pp.940-947
Minin, I. (Ed.). (March 2010). Microwave & Millimeter Wave Technologies Modern UWB
antennas and equipment, InTech, ISBN 978-953-7619-67-1, Croatia
Molisch A.F. , Foerster J. R. & Pendergrass, M. (2003). Channel models for ultrawide-band
personal area networks. IEEE Wireless Communications, Vol. 10, pp. 14–21
Palud S., Colombel F., Himdi M., Le Meins C. (2008). Compact multi-octave conical antenna.
Electronics Letters, Vol. 44, No. 11, pp 659-661
Sandler S. & King R.W.P. (1994). Compact conical antennas for wide-band coverage. IEEE
Transactions on Antennas & Propagation, Vol.42, No.3, pp 436-439
Sibille A., Roblin C., Bories S. & Lepage A. C. (2006). A Channel-Based Statistical Approach
to Antenna Performance in UWB Communications. IEEE Transactions on Antennas
& Propagation Vol. 54, No. 11, pp. 3207-3215
Stuzman W. A. & Thiele G. L. (1997). Antenna theory and Design, 2nd ed., John Wiley & Sons,
ISBN: 0-471-02590-9, New York
Taflove A. & Hagness S. C. (2005). Computational electrodynamics: the finite-difference time-
domain method, 3nd ed., pp. 411-472, ISBN 978-1-58053-832-9, Artech House, Boston
Yee K. (1996). Numerical solution of initial boundary value problems involving Maxwell's
equations in isotropic media. IEEE Transactions on Antennas and Propagation, Vol.14,
pp. 302–307
Yu Y. K. & Li J. (2008). Analysis of electrically small size conical antennas. Progress in
Electromagnetic Research Letters, Vol. 1, pp. 85-92
180 Scientific and Engineering Applications Using MATLAB
Zhou S., Ma J., Deng J., & Liu Q. (2009). A low-profile and broadband conical antenna.
Progress in Electromagnetic Research Letter, Vol.7, pp. 97-103
Wiesbeck W., Adamiuk G. & Sturm C. (2009). Basic properties and design principles of
UWB antennas. Proceedings of the IEEE, Vol. 97, No.2, pp.372-385
10
1. Introduction
A spherical shape of a submerged body with closed frame provides uniform drag at all
direction along its surface. In this chapter, the shape of a spherical URV that is used in this
book-chapter is presented. The vertical motion equation is also derived. The forces that
affect the dynamics of the system are also described. In order to control vertical motion due
to control depth position of the URV, a variable ballast mechanism is used. This mechanism
controls the weight of URV’s body. This chapter also presents detail mechanism of the
variable ballast system and describes detail of the used parts and design of the mechanism.
Kinematic and dynamic model of the variable ballast system are also derived.
2. Design of URV
The shape of spherical URV used in this book is shown in Figure 1. As a sphere body, the
location of center of buoyancy (COB) of URV’s body is at the center of sphere or the
intersection point between vertical and horizontal diameter. The variable ballast tank is
located at the top inside the hull. Location of the tank is adjusted so that the position of
center of mass (COM) is aligned vertically with COB. Mechanism of the variable ballast and
detail of its parts are explained in section 4.3. At the upper side of the hull above the tank,
there are some holes as the way of water to enter into and exit from the tank. The space
below the movable plate inside the hull is waterproofed so that the water can not enter this
space.
In order to make the URV stable in equilibrium condition, the hull of URV must be designed
with bottom heavy that is the center of mass is located at under of the equator or at
underside hemisphere. To make the hull in bottom heavy, fixed ballast is located at the
bottom of the hull. The location of COM of the hull must be aligned vertically with COB of
the hull thus in equilibrium condition, the position of the ballast tank is at the top of the
URV’s hull exactly. This condition is important when the URV is provided with horizontal
propulsion in order to give ability to the URV to move in horizontal plane.
182 Scientific and Engineering Applications Using MATLAB
W = FB + FD (1)
and,
W = mt g , (2.a)
FB = ρ wVfb g , (2.b)
1
FD = sign( v ) C D A fbρ w v 2 . (2.c)
2
The direction of gravitational force and buoyant force are opposite to each other when W is
downward and FB is upward. From Eq. 2.c, it can be seen that the direction of the drag force
depends upon the direction of the velocity. If the URV moves downward, the velocity is
positive so that the drag force is positive and it direction is upward. The drag force and
velocity are negative if the URV moves upward.
Variable Ballast Mechanism for Depth Positioning of a Spherical Underwater Robot Vehicle 183
1
mt g = ρ wVfb g + sign( v ) C D A fb ρ w v 2 . (3)
2
Since dimension of URV, Vfb and A fb , are constant, then the velocity v , is also constant.
This velocity is known as terminal velocity, which is expressed as
v = sign(mt g − ρ wVfb g )
(
2 mt g − ρ wVfb g ) , (4.a)
C D A fbρ w
v =
(
2 mt g − ρ wVfb g ) . (4.b)
C D A fbρ w
From Eq. 4.a, it can be seen that the vertical motion of the URV depends on the gravitational
force and the buoyant force. If W > FB , then the URV moves downward and it will move
upward if W < FB . If W = FB , the URV will stay at its position. Since volume of URV’s hull,
Vfb , g , and ρ w are constant, the buoyant force is also constant. So, the motion of URV
depends on the total mass of URV’s body, mt . By controlling mt , the vertical motion of the
URV can be controlled.
Since equilibrium condition is occurred when W = FB , then v = 0 and mt = ms which is
initial total mass of the URV. If the total mass changes as much as ∆m from the initial total
mass, ms , then the total mass of URV is expressed as
mt = ms + ∆m . (5)
Thus, the associated velocity will be change. The change of the velocity depends upon
whether ∆m is a variable or simply a constant. If ∆m is a variable, the acceleration, a ,
occurs. This acceleration, besides accelerates mass of URV itself, mt , also accelerates mass of
surrounding water which is known as added mass, ma .
184 Scientific and Engineering Applications Using MATLAB
f a = (ms + ∆m + m a )a . (6)
W = FB + FD + f a . (7)
1
ms g + ∆mg = FB + sign( v ) C D A fbρ w v 2 + (ms + ∆m + ma )a . (8)
2
Recalling equilibrium condition,
W = FB ,
mt = ms ,
v=0,
a=0. (9)
From Eq. 5 obviously we have ∆m = 0 .
By substituting ∆m and Eq. 9 into Eq. 8, yields
m s g = FB . (10)
1
∆m g = sign( v ) C D A fbρ w v 2 + (ms + ∆m + ma )a . (11)
2
Since ∆m g = ∆W , then Eq. 11 is written as
1 ∆W
∆W = sign( v ) C D A fbρ w v 2 + (ms + + ma ) a . (12)
2 g
By solving for the acceleration, a , the dynamic equation for vertical motion is given as (Xu
and Smith, 1994)
∆W sign( v )C D A fbρ w v 2
a= − . (13)
∆W ∆W
( ms + ma + ) 2(ms + ma + )
g g
And, if the depth position of the URV can be measured as z , then by differentiating z
respect to time t , the velocity of URV in vertical plane can be expressed as
v = zɺ . (14)
Variable Ballast Mechanism for Depth Positioning of a Spherical Underwater Robot Vehicle 185
Fig. 3.(a) Surface of water in ballast tank when the URV’s body is shaking; (b) Position
center of mass and center of buoyancy of water in the tank when tilt is change.
the movable plate is moving upward, the space of the tank will be decreased as well as the
volume of water in the ballast tank. If the movable plate is moving in opposite, downward,
the space of the tank will be increased and also the volume of water in the ballast tank.
Therefore, in any volume of water in the ballast tank there is no empty space in the ballast
tank that is not filled by water.
So, the screw converts the rotation motion into linear (vertical) motion. This coupling can be
seen in Figure 5(a). Based on Figure 5(a), l is lead of screw per revolution, ∆h denotes
change of nut position, ω3 is angular velocity of screw. If ∆t is the time needed by screw to
change nut position at ∆h regarding angular velocity ω3 , then their relation can be written
as
∆h l ω3
= . (15)
∆t 2 π
To turn the power screw, a DC motor is used and coupled with worm gear as illustrated in
Figure 5(b). The worm has number of thread per revolution equal to N w , and the gear has
number of teeth equal to N g . If the worm is coupled directly to the motor which turns in
velocity ωm , then the gear will turn in velocity ω2 which is expressed as
Nw
ω2 = ⋅ ωm . (16)
Ng
As shown in Figure 5(a), the gear and screw is ally so that its angular velocity is the same,
ω2 = ω3 . (17)
By substituting Eq. 17 and Eq. 16 into Eq. 15, the change of nut position can be rewritten as
∆h l N w
= ωm ,
∆t N g 2 π
l Nw
∆hɺ = ωm ,
N g 2π
N g 2 π ∆hɺ
ωm = . (18)
l Nw
T3 = TF + T fr , (19)
where TF is torque required to overcome force F , and Tfr is torque required to overcome
friction between screw and nut. To evaluate these terms, the equilibrium conditions are
applied such as illustrated in Figure 6.
188 Scientific and Engineering Applications Using MATLAB
Figure 6(a) illustrates coupling between nut and screw and also its parameters that must be
considered. There is an additional useful geometric relationship between lead angle, α , and
lead, l . Suppose the triangular segment of a plane wrapped around the screw is considered
in such a way that slanted edge lies along the helix and follows it for one revolution,
obviously we have
l
tan α = . (20)
π dm
Figure 6(b) illustrates a force P which is applied at a mean radius rm which causes the load
to be raised. The reactive forces act at point O on the screw thread surface.
The reactive force Fn acting normal to the surface has the following components:
OD = f r which is the friction force opposing movement up the thread surface
OA = is equal and opposite to the force being lifted. (F)
OB = is the vector sum of OD and OA and forms an angle θn with vector Fn
Summing the forces in the vertical direction results in
f r = µ s Fn . (22)
F
Fn = . (23)
cos θn cos α − µ s sin α
By equating Fn at Eq. 25 and Eq. 23, force P applied on screw in order to lift force F can be
expressed as
BC
tan θ n = = cos α tan θ . (27)
OB
Variable Ballast Mechanism for Depth Positioning of a Spherical Underwater Robot Vehicle 189
Fig. 6. (a) Screw and nut coupling; (b) Detail of forces working in the power screw
(roymech.co.uk, 2008)
If lead angle α is small, then cos α ≈ 1 , so we have
tan θn ≈ tan θ ,
θn ≈ θ . (28)
µ + cos θ tan α
P = F s . (29)
cos θ − µ s tan α
Again, by substituting Eq. 20 into Eq. 29, force P applied on screw to lift load F can be
rewritten as
190 Scientific and Engineering Applications Using MATLAB
πµ s dm + l cos θ
P = F . (30)
π dm cos θ − µ s l
In order to lift load F , a torque, T3U , must be applied to the screw. If the screw has mean
diameter dm , then the torque applied to the screw can be expressed as
dm
T3U = P . (31)
2
Substituting Eq. 30 into Eq. 31, the applied torque required to lift load F can be expressed as
dm πµ s dm + l cos θ
T3U = F . (32)
2 π dm cos θ − µ s l
In order to lower load F , a torque must be applied to the screw in reverse direction with
T3U and it is named as T3 L . Applying torque in reverse direction will also deliver force P in
reverse direction. By using same procedure in deriving T3U , the torque required to lower
load F can be expressed as
dm π ⋅ µ s ⋅ dm − l ⋅ cos θ
T3 L = F . (33)
2 π ⋅ dm cos θ + µ s ⋅ l
lw
tan λ w = . (34)
π dw
where Fwa is axial force of worm, Fwn is reactive force on worm, and Fs friction force of
worm. If coefficient friction of worm surface is µ w , then the friction force Fs is expressed as
Fs = µ w Fwn . (36)
If Eq. 36 is substituted into Eq. 35, then the reactive force Fwn is written as
Fwa
Fwn = . (37)
cos ϕn cos λ w − µ w sin λ w
If forces in horizontal direction are considered, then by summing of these forces will result
cos ϕn tan λ w + µ w
Fwt = Fwa . (39)
cos ϕn − µ w tan λ w
Then, by substituting Eq. 34 into Eq. 39, yields
192 Scientific and Engineering Applications Using MATLAB
l cos ϕn + πµ w dw
Fwt = Fwa w ,
π dw cos ϕn − µ w lw
π d cos ϕn − µ w lw
Fwa = Fwt w . (40)
lw cos ϕn + πµ w dw
From Figure 5(b), it is shown the relation of forces working at gear and worm. Forces on
gear are related by equilibrium to forces on the worm as
where Fgt and Fga are tangential and axial force working at gear respectively. If T2 is torque
applied on gear with diameter dg , then this torque T2 is expressed as
dg
T2 = Fgt . (42)
2
By equating Eq. 40 and Eq. 41.a and substitute into Eq. 42, yields
dg Fwt π dw cos ϕn − µ w lw
T2 = . (43)
2 lw cos ϕn + πµ w dw
To actuate this mechanism, the worm is coupled directly to the shaft of a DC motor. If Tm is
motor torque applied on worm to result tangential force Fwt which is expressed as
2
Fwt = Tm , (44)
dw
then by substituting Eq. 44 into Eq. 43, torque applied on gear is expressed as
dg Tm π dw cos ϕn − µ w lw
T2 = . (45)
dw lw cos ϕn + πµ w dw
Reviewing Figure 5 again, obviously can be seen that the gear and power screw are allied in
same shaft so that torque required to actuate the gear, T2 , will be equal to the torque
required to turn power screw, T3 . Since T3 = T2 , then T3U = T2U and T3U is torque needed by
screw to lift up the load F. In order to produce torque T3U on the screw or T2U on the gear,
the DC motor must produce torque TmU . If T3 = T3U , T2 = T2U , and Tm = TmU then by equating
Eq. 32 and Eq. 45 yields
dm dw F πµ s dm + l cos θ lw cos ϕn + πµ w dw
TmU = . (46)
2 d g π dm cos θ − µ s l π dw cos ϕn − µ w lw
By using the same analogy for calculating TmU , then the torque of the motor required to
lower the load F which is known as TmL , can be expressed as
dm dw F πµ s dm − l cos θ lw cos ϕn + πµ w dw
TmL = . (47)
2 dg π dm cos θ + µ s l π dw cos ϕn − µ w lw
From Eq. 46 and Eq. 47, it can be shown that many coefficients, which are constant, are
involved in the equation, so that if the constants are simplified then we have
dm dw
= kpr ,
2 dg
πµ s dm + l cos θ
= kTU ,
π dm cos θ − µ s l
πµ s dm − l cos θ
= kTL ,
π dm cos θ + µ s l
lw cos ϕn + πµ w dw
= kwg , (48)
π dw cos ϕn − µ w lw
where k pr is coefficient of power transmission ratio between worm gear set and power
screw, kTU and kTL are coefficient of power screw in lifting and lowering load mechanism
respectively, and kwg is coefficient of worm gear set. Hence, Eq. 46 and Eq. 47 can be
simplified into
or it can be written as
Tm = km F , (50)
where km = kmU = k pr kTU kwg ⇒ Tm = TmU and km = kmL = k pr kTL kwg ⇒ Tm = TmL .
From Eq. 50, it can be seen that torque Tm is the input, and F is the output. Although not
explicitly stated, it does not mean that if Tm = 0 then F must be zero. Because of friction, a
certain value of F must be reached to make it self-locking, before power screw start rotating
and allow the load lift or lower, and it’s called overhauling. To guarantee the screw will be
self-locking, a condition based on the geometric parameter and coefficient of friction must
be fulfilled (Collins, 2003).
194 Scientific and Engineering Applications Using MATLAB
If Figure 1 is analyzed, the force coming from inside hull, Fih , is caused by the change of air
pressure inside the hull, Pih , due to the change of space inside the hull. As explained before,
the variable ballast mechanism is used to control weight of URV by controlling volume of
water in ballast tank. To control volume of the water, a mechanism like piston is designed.
In this mechanism, a movable plate which is base of space in the tank that can be filled by
water is used. By controlling position of the movable plate which is height of the tank, the
volume of water in the tank can be controlled. Since space under movable plate is
impermeable, by the change of position of movable plate, the volume of space inside URV’s
hull is also change. This change impacts to the air pressure inside the hull.
As known that relation between pressure and volume, Vih , in closed space is expressed as
(Moran and Shapiro, 1998)
So, if the volume of air inside URV’s hull is changed, then its pressure is also changed. In
initial condition or in equilibrium condition, the pressure inside the hull is equal to the
pressure of the air at water surface, Pa . By assuming the air pressure at water surface and
Variable Ballast Mechanism for Depth Positioning of a Spherical Underwater Robot Vehicle 195
temperature inside URV’s hull are constant then since the volume of space inside the hull is
constant, the pressure inside the hull is also constant. If volume of the space inside the hull
is changed because of the change of position of movable plate, the pressure Pih is also
change and will cause a force act at movable plate surface, known as Fih . The relation of
Pa , Pih , and Fih is expressed as
where Avb is projected area of movable plate which is base of variable ballast tank.
In initial condition, where Pih = Pa , volume of the air or empty space inside URV’s hull is
Vih , and the position of movable plate is at the middle of full height of the tank, so if the
maximum height of the tank is h then the position of movable plate is at 0.5h from the top
of the tank, which is known as initial position. At this position, ∆h , which is the change of
movable plate position, is equal to zero ( ∆h = 0 ). So that if position of movable plate is
upper than initial position ( < 0.5h ) then ∆h < 0 and if lower than initial position then
∆h > 0 . By the change of ∆h , the volume of the air or space inside the hull will change as
∆V . The relation is expressed as
∆V = ∆h Avb . (53)
From Eq. 53, it is seen that the change of volume of the air inside the hull is equal to the
change of volume of water in the ballast tank.
By the change of volume of the air, the pressure is also change from its initial condition and
expressed as
Pa Vih
Pih = . (54)
Vih − ∆V
By substituting Eq. 53 into Eq. 54 and substitute the result into Eq. 52, then force acts on
movable plate’s surface due to change of volume of the air inside hull is expressed as
∆h ( Avb ) Pa
2
Fih = . (55)
Vih − ∆h Avb
From Figure 1, it can be shown that load F is resultant of Fhs , ∆W , and Fih , and the relation
is expressed as
where Wvb is weight of the water in the ballast tank and Fhs is hydrostatic force on surface
of movable plate. Weight of water in the ballast tank depends on volume of the water in this
tank. In equilibrium or initial condition, the volume of water is half of maximum volume of
the tank, Wvb = Wbs . By the change of position of movable plate in ∆h , the weight of water
in ballast is also will change in ∆W from initial weight. So that at any condition, Wvb can be
expressed as
196 Scientific and Engineering Applications Using MATLAB
Hydrostatic force, Fhs , is force acting on surface of immersed body caused by the height of
liquid (water) above of it, or in other word it can be said that hydrostatic force is weight of
the liquid above immersed surface. In this system, the height of liquid above is equal to the
depth position of the URV. If depth position of URV is measured form water surface to top
part of URV’s body known as z , then Fhs acting on surface of movable plate is expressed as
(Rajput, 2003)
∆h ( Avb ) Pa
2
∆W = ∆h Avb ρ w g ,
∆W
∆h Avb = . (60)
ρw g
Substituting Eq. 60 into Eq. 59, then the load F can be rewritten as
∆W Avb Pa
F = Wbs + ∆W + ρ w g Avb z −
ρ w gVih − ∆W . (61)
Recalling Eq. 50 and substitutes Eq. 61 into this equation, then torque of the motor that is
required to change position of movable plate in order to control amount of water in ballast
tank is expressed as
∆W Avb Pa
Tm = km Wbs + ∆W + ρ w g Avb z −
ρ w
gVih − ∆W . (62)
In order to change position of the movable plate, the DC motor must provide power Pm and
rotates at angular velocity ωm in order to produce torque at Tm , and can be expressed as
Pm = Tm ωm . (63)
∆W Avb Pa N g 2 π ∆hɺ
Pm = km Wbs + ∆W + ρ w g Avb z − . (64)
ρ w gVih − ∆W l N w
Variable Ballast Mechanism for Depth Positioning of a Spherical Underwater Robot Vehicle 197
N g 2π
If = k gc is known as coefficient of velocity reduction of power screw and worm gear
lNw
couple, then Eq. 64 can be rewritten as
∆W Avb Pa
Pm = k m k gc ∆hɺ W bs + ∆W + ρ w g Avb z − .
(65)
ρ w g Vih − ∆W
From Eq. 60, if A vb , ρ w , and g are simply constant, then by differentiating this equation
results
∆Wɺ
∆hɺ = , (66)
ρ w g Avb
where ∆hɺ is rate change of position of movable plate and ∆W ɺ is rate change of weight of
water in the ballast tank.
ɺ , then the rate change of weight of water
By substituting Eq. 66 into Eq. 65 and solving ∆W
in the ballast tank is obviously expressed as
ɺ = ρ w g A vb Pm . (67)
∆W
∆W A vb Pa
k m k gc W bs + ∆W + ρ w g A vb z −
ρ w g Vih − ∆W
saturation value depends on the maximum angular velocity of the DC motor that drives this
mechanism both in counterclockwise and clockwise direction. Then, the change of weight in
the ballast tank, ∆W , can be obtained by integrating ∆W ɺ which is shown in Figure 3. ∆W
also has saturation values ( ±∆Wsat ) which depends upon the maximum volume of the
ballast tank.
Output of this model is acceleration of the URV. In order to obtain velocity of URV’s vertical
motion, v , an integration block diagram is used. In order to get the depth position of the
URV, this velocity is integrated. The Simulink model is shown in Figure 5. The condition of
depth position and velocity are depth position always be positive ( z ≥ 0 ) and for
z=0⇒v≥0.
198 Scientific and Engineering Applications Using MATLAB
The Simulink model of URV motion in vertical plane which is taken from Eq. 13, is shown
in Figure 4. By combining all models, obviously Simulink model for depth positioning of the
spherical URV is shown in Figure 6.
Variable Ballast Mechanism for Depth Positioning of a Spherical Underwater Robot Vehicle 199
From Figure 6, it can be shown that the power input, Pm , has saturation values that is
± Pm max. This power depends upon the power provided by the DC motor.
URV’s hull:
Initial mass of URV ms 22.39 kg
Added mass of URV ma 11.2 kg
Diameter of URV D fb 0.35 m
Projected area of URV A fb 0.09616 m2
Initial volume of empty space inside URV’s hull Vih 50 % of Vfb
From Figure 7, it is also seen when the power is applied to the motor as positive value, the
weight change, ∆W , is increased and reaches saturation around 10N, as maximum value of
∆W . When the power is reset to zero, ∆W still remain at the last value, and the URV still
moves with velocity v , which is proportional to ∆W . The depth position of URV, z , will
increase since the velocity is available as positive.
If the velocity is negative, the URV is ascending as shown in Figure 8. The increment of
∆W depends on the total power given by the motor to actuate the variable ballast and also
depends on the depth position of the URV. If the power given is small then the increment of
∆W is also small, and since the power is available, ∆W will keep increasing till reach
saturation. If same power is given to the system but in different depth positions, the
increment of ∆W at deeper position is lower than shallower position. This is caused by the
availability of hydrostatic force.
When the power is zero, ∆W remains at its last value as well as velocity v . The velocity will
remain constant until ∆W change and velocity in this condition is known as terminal velocity.
This is the advantage of using variable ballast as vertical motion actuator, even the zero
power is given to the actuator, the URV still moves therefore it will save the power usage. If
∆W = 0 , then the zero velocity occurs. The depth position of URV, z , will remain at its last
position, and this condition is called equilibrium point. The equilibrium point occurs at any
depth position since v and ∆W is equal to zero.
If ramp input is applied, then the response of the system is shown in Figure 9. By looking to
the response, obviously the nonlinearity of the weight change of the URV’s body and the
velocity in vertical motion are shown. By the increment of Pm , v and ∆W also increase
until both of these reach saturation. Depth position, z , keeps increasing since v > 0 .
Variable Ballast Mechanism for Depth Positioning of a Spherical Underwater Robot Vehicle 203
6. Conclusion
The dynamic model of depth positioning of a spherical URV is obtained by considering the
physical laws involve in this system. The dynamic model presented in this chapter is
nonlinear. Some inputs are tested to see the responses and characteristics of this system. In
the next chapter, the control systems will be designed in order to control the depth position
of this spherical URV.
7. References
Collins, J. A. Mechanical Design of Machine Elements and Machines: a Failure Prevention
Perspective. New York: John Wiley & Sons Inc, 2003.
Moran, M. J., and H. N. Shapiro. Fundamentals of Engineering Thermodynamics, 3rd ed.
Chichester-England: John Wiley & Sons Inc, 1998.
204 Scientific and Engineering Applications Using MATLAB