System Modelling and Simulation
System Modelling and Simulation
System Modelling and Simulation
intentionally left
blank
Copyright © 2009, New Age International (P) Ltd., Publishers
Published by New Age International (P) Ltd., Publishers
Almost half a century has passed since System Analysis emerged as an independent field in Physical
Sciences. Number of books and research papers has appeared in the literature and a need is felt to
have a systematic one to the study of the subject. The basic techniques of Modeling and Simulation
are now being taught in undergraduate engineering courses, and its applications in various engineering
subjects require detailed studies. Similarly its application in “Weapon Systems Performance Analysis”
is very essential and has also been discussed in the book. Although no claim has been made about
the comprehensiveness of this book, yet attempt has been made to make this treatise useful to
engineers as well as scientists, especially defence scientists.
The present book is the output of my thirty years of work in the field of Armament and
System Analysis and also during my tenure in Institute of Engineering and Technology, Bhaddal,
where this subject is being taught at various levels. Most of the chapters in the book are based on
the papers published by the author in various technical journals. In order to make the analysis
easier to understand, basic mathematical techniques have also been discussed. It is not possible to
do full justice to the subject in one small book. But I have tried to condense as much material in
this book as possible. I will not say that this treatise is exhaustive, yet it may give quite insight
into the subject.
In the end I will like to acknowledge my friend and colleague Sh. Yuvraj Singh, Scientist E.,
Aeronautical Development Establishment, Bangalore, who was part of my team—Center for
Aeronautical System Studies and Analysis, Bangalore and had been quite helpful in preparing this
monograph.
V.P. Singh
This page
intentionally left
blank
CONTENTS
Preface v
0 . WHAT IS A SYSTEM 1– 7
1 . MODELING AND SIMULATION 9 – 25
1.1 PHYSICAL MODELS 10
1.2 MATHEMATICAL MODELS 12
1.2.1 Static Mathematical Models 13
1.2.2 Costing of a Combat Aircraft 13
1.2.3 A Static Marketing Model 15
1.2.4 Student Industrial Training Performance Model 16
1.3 COMPUTER MODELS 18
1.3.1 Runway Denial using BCES Type Warhead 18
1.3.2 Distributed Lag Models—Dynamic Models 19
1.4 COBWEB MODELS 20
1.5 SIMULATION 23
1.5.1 Monte Carlo Simulation 24
2.3.2 Variance 35
2.3.3 Some Theorems on Variance 35
2.4 MEASURE OF PROBABILITY FUNCTION 36
2.4.1 Central Tendency 36
2.4.2 Median 36
2.5 SOME IMPORTANT DISTRIBUTION FUNCTIONS 37
2.5.1 Cumulative Distribution Function 37
2.5.2 Uniform Distribution Function 38
2.5.3 Binomial Distribution Function 39
2.5.4 Poisson’s Distribution 44
2.6 CONTINUOUS RANDOM VARIABLE 46
2.7 EXPONENTIAL DISTRIBUTION 48
2.7.1 Gamma Distribution Function 49
2.7.2 Erlang Density Function 50
2.8 MEAN AND VARIANCE OF CONTINUOUS DISTRIBUTION 51
2.9 NORMAL DISTRIBUTION 52
2.9.1 Properties of Normal Distribution Curve 53
2.9.2 Cumulative Density Distribution Function of Normal Distribution 53
2.9.3 An Experiment for the Demonstration of Normal Distribution Function 55
2.9.4 Example of Dispersion Patterns 56
2.9.5 Estimation of Dispersion 56
2.10 CIRCULAR PROBABLE ERROR (CEP) AND THE PROBABLE ERROR (PE) 58
2.10.1 Range and Deflection Probable Errors 59
2.10.2 Probability of Hitting a Circular Target 60
APPENDIX 2.1 63
BIBLIOGRAPHY 241
INDEX 245
This page
intentionally left
blank
123456789
123456789
123456789
123456789
123456789
0
123456789
123456789
123456789
WHAT IS A SYSTEM
System modeling and computer simulation, recently has become one of the premier subject in the
industry as well as Defence. It helps an engineer or a scientist to study a system with the help of
mathematical models and computers. In other words, one can say, system simulation is nothing but
an experiment with the help of computers, without performing actual experiment. It saves lot of money
which is required, if we actually perform experiments with the real system.
Present book, “System Modeling and Simulation” has been written by keeping engineering students
as well as scientists especially defence scientists in mind. Earlier this manuscript was prepared only for
the use of defence scientists which has now been extended to other engineering applications. Modeling
of a weapon system is just an application of basic techniques of modeling and simulation, which are
being discussed in chapter two and four. After my superannuation from Defence Research & Development
Organisation, briefly called DRDO in 2000, when I joined the Punjab Technical University, and taught
this subject to B. Tech and M. Tech students, the manuscript was rewritten, so that it should be useful
to engineering students too. Although many of the examples have been taken from target damage, yet
care has been taken to include other examples, from marketing, and mechanical engineering and other
related subjects. My intentions are that this book should prove itself as a complete manual for system
simulation and modeling. That is the reason, that basic subjects required for modeling and simulation,
such as probability theory, numerical methods and C++ have also been included and discussed wherever
required. Wherever possible, computer programmes have been given with the output.
First question the user of this book can raise is, after all what is a system1? We give here a
popular definition of what a system is? A broader definition of a system is, “Any object which has
some action to perform and is dependent on number of objects called entities, is a system”. For
example a class room, a college, or a university is a system. University consists of number of colleges
(which are entities of the system called university) and a college has class rooms, students,
laboratories and lot many other objects, as entities. Each entity has its own attributes or properties.
For example attribute of a student is to study and work hard. Each college in itself can be treated
as a complete system. If we combine few of these objects, joined in some regular interactions or
inter-dependence, then this becomes a large system. Thus we can say university is a large system
whereas college is a system. This means, each system can be divided into blocks, where each block
is in itself a complete and independently working system (Fig. 0.1). When these blocks are combined,
depending on some interdependence, they become entities of a larger system. An aircraft for example,
is another example of a system. It consists of a cockpit, pilot, airframes, control system, fuel tank,
engines etc. Each of these blocks can in itself be treated as a system and aircraft consisting of
these interdependent blocks is a large system. Table 0.1 gives components of a system (say college)
for illustrations purpose.
Does a system also has some attributes? I will say yes. Systems broadly can be divided into two
types, static system and dynamic system. If a system does not change with time, it is called a Static
System and if changes with time, it is called a Dynamic System. Study of such a system is called
System Analysis. How this word originated is of interest to know.
COLLEGE
PLAY CONFERENCE
RAW STUDENTS DEPARTMENTS CANTEEN HOSTELS GRADUATE
GROUND HALL
While looking at these systems, we see that there are certain distinct objects, each of which
possesses some properties of interest. There are also certain interactions occurring in the system
that cause changes in the system. A term entity will be used to denote an object of interest in a
system and the term attributes denotes its properties. A function to be performed by the entity is
called its activity. For example, if system is a class in a school, then students are entities, books
are their attributes and to study is their activity. In case of the autopilot aircraft discussed below,
entities of the system are gyroscope, airframe and control surfaces. Attributes respectively are
gyroscope setting, speed and control surface angles. Activity of the aircraft is to fly. Banks et al.,
(2002) has defined state variables and events as components of a system. The state is defined as
collection of variables necessary to describe the system at any time, relative to the objectives of
the study. An event is defined as an instantaneous occurrences that may change the state of the
system. But it is felt, entities are nothing but state variables and activity and event are similar. Thus
there is no need of further bifurcation.
Sometimes the system is effected by the environment. Such a system is called exogenous. If it
is not effected by the environment, it is called endogenous. For example, the economic model of a
country is effected by the world economic conditions, and is exogenous model. Aircraft flight is
exogenous, as flight profile is effected by the weather conditions, but static model of the aircraft is
endogenous. A class room in the absence of students, is endogenous. As mentioned earlier study of
a system is called System Analysis. How this word “System Analysis” has cropped up? There is an
interesting history behind this.
What is a System 3
After second world war (1939–1945), it was decided that there should be a systematic study
of weapon systems. In such studies, topics like weapon delivery (to drop a weapon on enemy) and
weapon assignment (which weapon should be dropped?) should also be included. This was the time
when a new field of science emerged and the name given to it was “Operations Research”. Perhaps
this name was after a British Defence Project entitled “Army Operations Research” [19]*. Operations
Research at that time was a new subject. Various definitions of “Operations Research” have been given.
According to Morse and Kimball [1954], it is defined as scientific method of providing executive
departments with a quantitative basis for decisions regarding the operations under one’s control.
According to Rhaman and Zaheer (1952) there are three levels of research : basic research, applied
research and operations research. Before the advent of computers, there was not much effect of this
subject on the weapon performance studies as well as other engineering studies. The reason being
that enormous calculations were required for these studies, which if not impossible, were quite difficult.
With the advent of computers, “Computerised Operations Research” has become an important subject
of science. I have assumed that readers of this book are conversant with numerical methods as well
as some computer language.
History of Operations Research has beautifully been narrated by Treften (1954). With the time,
newer and newer techniques were evolved and name of Operations Research also slowly changed
to “System Analysis” [17]. Few authors have a different view that “Operations Research” and “System
Analysis” are quite different fields. System Analysis started much later in 1958–1962 during the
Kennedy administration in US (Treften, FB 1954). For example, to understand any system, a scientist
has to understand the physics of the system and study all the related sub-systems of that system.
In other words, to study the performance evaluation of any system, one has to make use of all the
scientific techniques, along with those of operations research. For a system analyst, it is necessary
to have sufficient knowledge of the every aspect of the system to be analysed. In the analysis of
performance evaluation of a typical weapon as well as any other machine, apart from full knowledge
of the working of the system, basic mathematics, probability theory as well as computer simulation
techniques are used. Not only these, but basic scientific techniques, are also needed to study a system.
In the present book our main emphasis will be on the study of various techniques required for
system analysis. These techniques will be applied to various case studies with special reference
to weapon system analysis and engineering. Modeling and simulation is an essential part of system
analysis. But to make this book comprehensive, wherever possible, other examples will also be given.
Present book is the collection of various lectures, which author has delivered in various courses
being conducted from time to time for service officers and defence scientists as well as B. Tech
* Number given in the square bracket [1], are the numbers of references at the end of this book. References however are
in general given as author name (year of publication).
4 System Modeling and Simulation
and M. Tech students of engineering. Various models from other branches of engineering have also
been added while author was teaching at Institute of Engineering and Technology, Bhaddal, Ropar
in Punjab. Attempt has been made to give the basic concepts as far as possible, to make the treatise
easy to understand.
It is assumed that student, reading this book is conversant with basic techniques of numerical
computations and programming in C++ or C language. But still wherever possible, introductory sections
have been included in the book to make it comprehensive.
As given earlier, Simulation is actually nothing but conducting trials on computer. By simulation
one can, not only predict the system behaviour but also can suggest improvements in it. Taking any
weapons as an example of a system, lethal capabilities of the weapon which are studied by evaluating
its terminal effects i.e., the damage caused by a typical weapon on a relevant target when it is dropped
on it, is also one of the major factor. There are two ways of assessing the lethal capabilities of the
weapons. One is to actually drop the weapon on the target and see its end effects and the other is
to make a model using a computer and conduct a computer trial, which simulates the end effects
of the weapon on the target. Conducting trials on the computers is called computer simulation. Former
way of conducting trial (actual trials) is generally a very costly affair. In this technique, cost of the
weapon as well as target, are involved. But the second option is more cost effective.
As an example of a conceptually simple system (Geofrey Gordon, 2004), consider an aircraft
flying under the control of an autopilot (Fig. 0.2). A gyroscope in the autopilot detects the difference
between the actual heading of aircraft and the desired heading. It sends a signal to move the control
surfaces. In response to control surfaces movement, the aircraft steers towards the desired heading.
θi ε CONTROL
GYROSCOPE AIRFRAME
DESIRED SURFACES
HEADINGS
θο ACTUAL HEADING
A factory consisting of various units such as procurement department, fabrication and sale
department is also a system. Each department of the factory, on which it depends, is an independent
system and can be modelled independently. We have used the word “modelled” here. What is modeling,
will be discussed in chapter one of this book.
Scientific techniques used in system studies can also broadly be divided into two types:
1. Deterministic studies and
2. Probabilistic studies
Deterministic studies are the techniques, where results are known exactly. One can represent
system in the form of mathematical equations and these equations may be solved by analytic methods
or by numerical computations. Numerical computations is one of the important tools in system analysis
and comes to rescue the system analyst, where analytical solutions are not feasible. In case of study
of damage to targets, knowledge of shock waves, fragments penetration, hollow charge and other
allied studies is also required. Some of the topics have been introduced in this book for the purpose
of broader horizon of the student’s knowledge.
What is a System 5
Apart from the two types of studies given above, system can be defined as
(i) Continuous and
(ii) Discrete.
Fluid flow in a pipe, motion of an aircraft or trajectory of a projectile, are examples of continuous
systems. To understand continuity, students are advised to refer some basic book on continuity.
Examples of discrete systems are, a factory where products are produced and marketed in lots. Motion
of an aircraft is continuous but if there is a sudden change in aircraft’s level to weather conditions,
is a discrete system. Another example of discrete system is firing of a gun on an enemy target.
It is important to conduct experiments to confirm theoretically developed mathematical models.
Not only the experiments are required for the validation of the theoretical models but various parameters
required for the development of theoretical models are also generated by experimental techniques. For
example to study the performance of an aircraft, various parameters like drag, lift and moment
coefficients are needed, which can only be determined experimentally in the wind tunnel. Thus theory
and experiment both are complementary to each other and are required for correct modeling of a
system. In case of marketing and biological models in place of experiments, observations and trends
of the system over a time period, are required to be known for modeling the system. This issue will
be discussed in chapter eight.
The aim of this book is to discuss the methods of system analysis vis a vis an application to
engineering and defence oriented problems. It is appropriate here to make reference to the books, by
Billy E. Gillett (1979), Pratap K. Mohapatra et al., (1994) and Shannon (1974), which have been
consulted by the author.
Layout of the book is arranged in such a way that it is easy for the student to concentrate on
examples of his own field and interest. Our attempt is that this book should be useful to students of
engineering as well as scientists working in different fields. Layout of the book is as follows:
In chapter one, basic concepts of modeling and simulation are given. Number of examples are
given to illustrate the concept of modeling. Also different types of models are studied in details.
In chapter two, basic probability theory, required for this book has been discussed. Probability
distribution functions, as used in modeling of various problems have been included in this chapter.
Wherever needed, proof for various derivations are also included. Application of probability theory
to simple cases is demonstrated as examples. Although it is not possible to include whole probability
theory in this book, attempts have been made to make this treatise comprehensive.
Chapter three gives a simple modeling of aircraft survivability, by overlapping of areas, using
probability concepts developed in chapter two. It is understood that projection of various parts on
a plane are known. This chapter will be of use to scientists who want to learn aircraft modeling tools.
However engineering students, if uncomfortable, can skip this chapter.
Simulation by Monte Carlo technique using random numbers, is one of the most versatile
technique for discrete system studies and has been dealt in chapter four. The problems which
cannot be handled or, are quite difficult to handle by theoretical methods, can easily be solved
by this technique. Various methods for the generation of uniform random numbers have been
discussed. Properties of uniform random numbers and various tests for testing the uniformity of
random numbers are also given. Normal random numbers are of importance in various problems
in nature, including weapon systems. Methods for generating normal random numbers have been
discussed in details. Study of different types of problems by computer simulation technique have
been discussed in this chapter. Here simulation technique, and generation of different types of
6 System Modeling and Simulation
random numbers is studied in details. Computer programs in C++ for different techniques of
random number generation are also given.
Chapter five deals with the simulation and modeling of Continuous Dynamic Systems. Problem
becomes slightly complex, when target is mobile. Problem of mobile targets and surface to air weapons,
will be studied in chapter six. A study of vulnerability of aerial targets such as aircraft will be discussed.
Simulation and modeling can not be learnt without working on practical problems. Therefore number
of examples have been worked out in this chapter.
A case study of aircraft survivability has been given in chapter six. Aircraft model discussed
in chapter three was purely on probability concepts. But present model is a simulation model, where
all the techniques developed in first five chapters has been utilised. What is the probability that an
aircraft will survive when subjected to ground fire in an enemy area, has been studied in this model.
This model involves mathematical, and computer modeling. Concept of continuous and stochastic
modeling has also been used wherever required.
Simulation of manufacturing process and material handling is an important subject of Mechanical
Engineering. Queuing theory has direct application in maintenance and production problems. Queuing
theory has been discussed in chapter seven. Various applications on the field of manufacturing process
and material handling have been included in this chapter.
There are various phenomena in nature which become unstable, if not controlled in time. This
is possible only if their dynamics is studied and timely measures are taken. Population problem is one
of the examples. Nuclear reaction is another example. This type of study is called Industrial dynamics.
Eighth chapter deals with System Dynamics of such phenomena. Various cases of growth and decay
models with examples have been discussed in this chapter.
One of the pressing problems in the manufacturing and sale of goods is the control of inventory
holding. Many companies fail each year due to the lack of adequate control of inventory. Chapter nine
is dedicated to inventory control problems. Attempt has been made to model various inventory control
conditions mathematically.
Costing of a system is also one of the major job in system analysis. How much cost is involved
in design and manufacturing of a system is discussed in tenth chapter. Cost effectiveness study is very
important whether procuring or developing any equipment. In this chapter basic concepts of costing
vis a vis an application to aircraft industry has been discussed.
EXERCISE
This page
intentionally left
blank
Modeling and Simulation 9
12345678
12345678
12345678
12345678
12345678
1
12345678
12345678
12345678
Term Modeling and Simulation is not very old and has become popular in recent years. But the field of
modeling and simulation is not new. Scientific workers have always been building models for different
physical scenarios, mathematically as well as in laboratories. Here question arises, what is the basic
difference between modeling and simulation? In fact, there is very thin border line between the both. In
fact, physical models can not be called simulation, but mathematical models can be called simulation.
Similarly computer simulation can not be given name of modeling but mathematical simulation can be
called modeling. Model of a system is the replica of the system, physical or mathematical, which has all
the properties and functions of the system, whereas simulation is the process which simulates in the
laboratory or on the computer, the actual scenario as close to the system as possible. In fact, a modeling
is the general name whereas simulation is specific name given to the computer modeling. Models can be
put under three categories, physical models, mathematical models and computer models. All of these
types are further defined as static and dynamic models. Physical model is a scaled down model of actual
system, which has all the properties of the system, or at least it is as close to the actual system as
possible. Now-a-days small models of cars, helicopters and aircraft are available in the market. These
toys resemble actual cars and aircraft. They are static physical models of dynamic systems (cars,
helicopters and aircraft are dynamic systems). In wind tunnel, scaled down models of an aircraft are
used to study the effect of various aerodynamic forces on it. Similarly before the construction of big
buildings, scaled down models of the buildings are made. Well known laws of similitude are used to
make the laboratory models. All the dimensions are reduced with respect to some critical lengths
(Sedov LI, 1982). Figure 1.1 gives different types of models. These types will be studied in details in the
coming chapters.
While building a model certain basic principles are to be followed. While making a model one
should keep in mind five basic steps.
• Block building
• Relevance
• Accuracy
• Aggregation
• Validation
10 System Modeling and Simulation
A model of a system can be divided into number of blocks, which in itself are complete systems.
But these blocks should have some relevance to main system. For example, let us take an example of
a school. Class rooms are blocks of the school. Aim of the school is to impart education to students
and class rooms are required for the coaching. Thus relevance of class rooms (blocks) with school
is coaching. Interdependency is the one of the important factor of different blocks. Each block should
be accurate and tested independently. Then these blocks are to be integrated together. Last is the validation
i.e., this model is to be tested for its performance. For validation following methods can be used.
• If the model is mathematical model, then some trivial results can be run for verifications.
• If experimental results are available, model can be checked with these experimental results.
In the following sections, we will discuss in details the various types of models as shown in Fig. 1.1.
MODEL
SYSTEM SIMULATION
Let us take an example of hanging wheel of a stationary truck and analyze its motion under various
forces. Consider a wheel of mass M, suspended in vertical direction, a force F(t), which varies with
time, is acting on it. Mass is connected with a spring of stiffness K, and a piston with damping factor
D. When force F (t), is applied, mass M oscillates under the action of these three forces. This model
can be used to study the oscillations in a motor wheel. Figure 1.2 shows such a system. This is
a discrete physical static model. Discrete in a sense, that one can give discrete values F and observe
the oscillations of wheel with some measuring equipment. When force is applied on it, which is a
function of time, this discrete physical static model becomes dynamic model. Parameters K and D
can also be adjusted in order to get controlled oscillations of the wheel. This type of system is called
spring-mass system. Load on the beams of a building can be studied by the combination of spring-
mass system. Mathematical model of this system will be studied in coming chapters.
Let us consider another static physical model which represents an electric circuit with an
inductance L, a resistance R, and a capacitance C, connected with a voltage source which varies
with time, denoted by the function E (t). This model is meant for the study of rate of flow of current
as E (t) varies with time. There is some similarity between this model and model of hanging wheel.
It will be shown below that mathematical model for both is similar. These physical models can easily
be translated into a mathematical model. Let us construct mathematical model of system describing
hanging wheel. Using Newton’s second law of motions, system for wheel model can be expressed
in the mathematical form as
d2 x dx
M 2 + D + Kx = KF(t) ...(1.1)
dt dt
where x = the distance moved,
M = mass of the wheel,
K = stiffness of the spring,
D = damping force of the shock absorber.
This system is a function of time and is a dynamic system. Equation (1.1) cannot be solved
analytically and computational techniques can be used for solving the equations. Once we interpret
the results of this equation, this becomes a dynamic mathematical model of the system. Physical model
12 System Modeling and Simulation
can also be used to study the oscillations by applying force F, and measuring the displacement x of
the model. Values of K and D can be changed for required motion with minimum oscillations. This
equation will be discussed in the section on Mathematical models in section 1.2. Equation of electrical
circuit given above can be written as
q E (t )
Lq + Rq + =
C C
.
Here the q is the electric current. In this equation if we say q = x, L = M, R = D, E = F and
K = 1/C, then one can easily see that this equation is similar to (1.1).
R L
E(t) C
Thus same mathematical model, by using different constants can give solution for hanging wheel
as well as electrical circuit.
d2 x dx
+ ω 2 x = ω F (t )
2
+ 2ζω ...(1.2)
dt 2 dt
where 2ζω = D / M and ω 2 = K /M. Expressed in this form, solution can be given in terms of the
variable ωt. Details of the solution of this equation will be given in chapter five on continuous models.
Expressed in this form, solution can be given in terms of the variable ωt , where ω is the frequency
of oscillation given by
ω2 = 2π f
and f is the number of cycles per second.
Modeling and Simulation 13
We can integrate equation (1.2) using numerical techniques. It will be shown that wheel does
not oscillate for ζ ≥ 1 .
the target and due to failure of weapon release mechanism respectively. Then Na, Nb and Nc are given
as follows:
Na = N . pa
Nb = N . (1 – pa ) pb ...(1.4)
Nc = N . (1 – pa ) (1 – pb ) pc
In equation (1.4), Nb is the product of aircraft not aborted due to malfunctioning i.e., N (1 – pa)
and probability of abort due to missing the target ( pb ). Same way Nc is evaluated. Since Na number
of aircraft have not flown for the mission, we assume, that they are not subjected to enemy attrition.
Then total number of aircraft going for mission is
N1 = N – Na = N (1 – pa ) ...(1.5)
Nb and Nc are number of aircraft being aborted after entering into enemy territory and thus are
subject to attrition. Therefore total number of aircraft accomplishing the successful mission is
N – Na – Nb – Nc
Out of which surviving aircraft are given as
Ns = p(N – Na – Nb – Nc )
or by using (1.4) in the above equation,
Ns = p .N [1 – pa – (1 – pa ) pb – (1 – pa ) (1 – pb ) pc ] ...(1.6)
There are another category of surviving aircraft (given by equation (1.4)) which have not been
able to accomplish the mission. Amongst these Na aircraft aborted before entering the enemy territory
and are not subject to attrition, where as Nb and Nc category have attrition probability p, as they abort
after reaching near the target. Let these aircraft be denoted by Nas. Then
Nas = Na + ( Nb + Nc ) . p ...(1.7)
Therefore total number of survived aircraft is
Ns + Nas
Once number of survived aircraft is evaluated, we can fix a cost factor. This is accomplished
as follows.
If the acquisition cost of an aircraft, whose total life in H hours is C, then the life cycle cost
of the aircraft, is given by
C (1 + k)
where k is the cost factor for computation of life cycle cost which includes cost of spares and
maintenance, operation etc., and strike-off wastage during its life. Value of cost for computation of
life cycle cost which includes cost of various aircraft is given in Appendix 10.1. The cost of the
surviving aircraft for the assigned mission, and have performed the mission is given by
C1 = (C (1 + k) Ns t/H ), ...(1.8)
where t is mission sortie time in hours.
where ta, tb, tc are the sortie time taken in each of the aborted cases mentioned above. Timings ta, tb
and tc are input data and is based on the war estimates and total sortie time t of the aircraft. For example
tc = t as in this case the aircraft had reached the target but could not release the weapon due to failure
of firing mechanism. Similarly tb can also be taken as equal to t as aircraft has flown up to the target
but could not locate it.
(c ) Cost of killed aircraft
Let Nk be the total number of killed aircraft given by Nk = (N – Na) . (1 – p) and tj is the life already
spent by Nj-th aircraft, then the total cost of killed aircraft is given by
p
C Csow
C3 = 1 + ∑ Nj ( H − tj ) ...(1.10)
H C j =1
where Nk = N1 + N2+ …+ Np and the terms Csow / C in above equation has been taken due to the
reasons that when an aircraft which has spent a life of tj hours is lost, partial cost of spares, fuel,
and maintenance is not lost. Only loss in this case is C + Csow. Thus if we assume that all the aircraft
are new i.e., tj = 0, we have,
Csow
C3 = N k C 1 + ...(1.11)
C
It is to be noted that in this equation only cost due to strike-off wastage is taken out of cost factor
k. This is because cost of spares and maintenance is not assumed as gone when an aircraft is killed.
which demand can match the supply. Let us model this situation mathematically. If we denote price
by P, supply by S and demand by D, and assuming the price equation to be linear we have
D = a – bP
S = c + dP ...(1.15)
S = D
In the above equations, a, b, c, d are parameters computed based on previous market data.
Equation S = D says supply should be equal to demand so that market price should accordingly
be adjusted. Let us take values of a = 500, b = 2000, c = –50 and d =1500. Value of c is taken
negative, since supply cannot be possible if price of the item is zero. In this case no doubt equilibrium
market price will be
a−c 550
P = = = 0.1571
b+d 3500
and S = 186
Price
0.25
Demand Supply
Price
0.15
0.10
0.05
In this model we have taken a simplistic linear case but equation (1.15) may be complex. In that
case solution may not be so simple. More usually, the demand and supply are depicted by curves
with slopes downward and upward respectively (Fig. 1.5). It may not be possible to express the
relationships by equations that can be solved easily. Some numerical or graphical methods are used
to solve such relations. In addition, it is difficult to get the values of the coefficients of the model.
Observations over the extended period of time, however, will establish the slopes (that is values of
b and d ) in the neighbourhood of the equilibrium points. These values will often fluctuate under the
global and local economic conditions.
1.2.4 Student Industrial Training Performance Model
For engineering students, six months training in industry is a part of their course curriculum and is
compulsory. It has been observed that training marks allotted to students from industrial Institutions,
vary drastically, irrespective of the academic record of the student. Nature of project offered and
standard of training institutions pay a very dominating role in these criteria. Due to this sometimes,
very good students, who are supposed to top the University exam, can suffer with no fault on their
part. A model to optimize the industrial marks which are about 40% of total marks is presented below.
Let M be the marks allotted by industry and M the optimized marks.
Modeling and Simulation 17
Then,
M = 0.75 + 0.5 × (0.15) + 0.25 × (0.10) + 0.125 × (0.15) – 0.125 × (0.6 – 0.8)
⇒ M = 0.75 + 0.075 + 0.025 + 0.019 + 0.025
∴ M = 0.894
(II) Case Study (Poor Student)
M = 0.95
mac = 0.50
mext = 0.6
mint = 0.5
m ind = 0.9
18 System Modeling and Simulation
Then,
M = 0.95 + 0.5 × (0.5 – 0.8) + 0.25 × (0.6 – 0.8) + 0.125 × (0.5 – 0.8) – 0.125 × (0.9 – 0.8)
⇒ M = 0.95 + 0.15 – 0.05 – 0.034 – 0.0125
∴ M = 0.70
Model can be used for upgrading or down grading (if required) the industrial marks.
SAMPLE DMAI
Fig. 1.6: A typical airfield, showing runways and DMPIs (black sectors).
To achieve this, we draw certain areas on the track so that if atleast one bomb falls on each
of these areas, 1000 meters will not be available anywhere on the three tracks. These areas are called
Desired Mean Areas of Impact (DMAI ), and have been shown as black areas in Fig. 1.6. Desired
Mean Points of Impact (DMPIs) and strips are chosen in such a way that, if each strip has at least
one bomb, no where a strip of dimensions 1000m × 15m is available. Number of strips Ns of effective
width Ws in a DMAI is given by,
1, if Wd = W
Ns = 2W ...(1.16)
int W + 2r + 1, otherwise
d b
where W, Wd , rb are the width, denial width respectively of the RW and lethal radius of the bomblet.
Monte Carlo computer model of airfield denial has been discussed by Singh et al., (1997) and
is given in chapter four.
Then
C = 20 + 0.7(Y − T )
I = 2 + 0.1Y
T = 0 + 0.2Y ...(1.17)
Y =C+ I +G
All quantities are expressed in billions of rupees.
This is a static model, but it can be made dynamic by picking a fixed time interval, say one year,
and expressing the current values of the variables in terms of values of the previous year. Any variable
that appears in the form of its current value and one or more previous year’s values is called lagged
variables. Value of the previous year is denoted by the suffix with-1.
The static model can be made dynamic by lagging all the variables, as follows;
C = 20 + 0.7(Y−1 − T−1 )
I = 2 + 0.1Y–1
T = 0.2Y−1 ...(1.18)
Y = C−1 + I−1 + G−1
In these equations if values for the previous year (with –1 subscript) is known then values for
the current event can be computed. Taking these values as the input, values for the next year can
also be computed. In equation (1.18) we have four equations in five unknown variables.
It is however not necessary to lag all the variable like it is done in equation (1.18). Only one
of the variable can be lagged and others can be expressed in terms of this variable. We solve equation
for Y in equation (1.17) as
Y = 20 + 0.7(Y – 0.2Y ) + I + G
= 20 + 0.56Y + I + G
or Y = 45.45 + 2.27(I + G)
Thus we have,
I = 2.0 + 0.1Y–1
Y = 45.45 + 2.27( I + G )
T = 0.2Y ...(1.19)
C = 20 + 0.7(Y − T )
In equations (1.19) only lagged parameter is Y. Assuming that government expenditure is known
for the current year, we first compute I. Knowing I and G, Y and T for the current year is known,
and thus C is computed from the last equation. In this problem, lagged model is quite simple and
can be computed with hand calculator. But national economic models are generally not that simple
and require long computations with number of parameters.
price and demand of a product in the market subject to a condition that supply and demand should
be equal. But supply of the product in the market depends on the previous year price, and that can
be taken as lagged parameter. Thus equations (1.15) become
D = a – bP
S = c + dP–1 ...(1.20)
D = S
Model 1 Model 2
P 0 = 30 P0 = 5
a = 12 a = 10.0
b = 30 b = 0.9
c = 1.0 c = –2.4
0.9 1.2
In the equations (1.20), values of parameters a, b, c, and d can be obtained by fitting a linear
curve in the data from the past history of the product. We assign some initial value to the product
say P0 and find S from second equation of (1.20). Thus S and D are known and first equation of
(1.20) gives us new value of P. Using this value of P as initial value, we repeat the calculations and
again compute P for the next period. If price converges, we say model (1.20) is stable. Let us take
two examples and test whether these models converge or not.
Model 1 Model 2
i P i P
0 –0.533333 0 7.11111
1 0.382667 1 4.2963
2 0.355187 2 8.04938
3 0.356011 3 3.04527
4 0.355986 4 9.71742
5 0.355987 5 0.821216
6 0.355987 6 12.6828
7 0.355987 7 –3.13265
8 0.355987 8 17.9546
9 0.355987 9 –10.1618
10 0.355987 10 27.3268
We have given a small program to find the value of P on the next page.
22 System Modeling and Simulation
10 Demand
9
8
7 1
6 3 5 Supply
Price
5
6
4 4
3 2
2
1 0
0
1 2 3 4 5 6 7 8 9 10
Quantity
For simulating the real life systems, no real science can serve the purpose, because knowledge
of different discipline is generally needed. That is why sometimes simulation is called an art and not
science. In the coming chapters, we will study different techniques, required for simulation. Probability
theory is one of the important scientific fields required for stochastic simulation. In next chapter, we
will study probability in details.
1.5 SIMULATION
In the section 1.2, we had given an example of a mathematical model of hanging wheel of a vehicle.
It will be shown by numerical computations of the equation (1.1) in chapter five, that system does
not oscillate when parameter ζ is greater than or equal to one. This we could find by numerically
integrating the equation (1.2). If it was possible to get analytical solution of equation (1.2), one could
easily find by putting ζ = 1, that system does not oscillate. However, with the method of numerical
techniques, one has to run the program for different values of ζ and then find out that for ζ ≥ 1 , system
is stable. Same is the case with simulation. One has to run simulation model number of time, to arrive
at a correct output. Due to this reason, sometimes numerical computation is called simulation. But
information derived numerically does not constitute simulation. Numerical computation only gives
variations of one parameter in terms of other and does not give complete scenario with the time.
Simulation has long been an important tool of designers, whether they are simulating a supersonic
jet, a telephone communication system, a wind tunnel, a large scale military war gaming, or a
maintenance operation.
Although simulation is often viewed as a “method of last resort” to be employed when every other
technique has failed. Recent advances in simulation methodologies, availability of softwares, and
technical developments have made simulation one of the most widely used and accepted tools in system
analysis and operation research.
Naylor et al., [41] defines the simulation as follows:
Simulation is a numerical technique for conducting experiments on a digital computer, which
involves certain types of mathematical and logical models over extended period of real time.
We thus define system simulation as the technique of solving problems by the observation of the
performance, over time, of a dynamic model of the system. In other words, we can define simulation
as an experiment of physical scenario on the computer.
Why simulation is required? According to Naylor [41], some of the reasons why simulation is
appropriate are:
1. Simulation makes it possible to study and experiment with the complex internal interactions
of a given system, whether it be a firm, an industry, an economy, or some subsystem
of one of these.
2. Through simulation we can study the effect of certain informational, organizational, and
environmental change on the operation of a system by making alterations in the model of
the system and observing the effects of these alterations on the system’s behavior.
3. Detailed observation of the system being simulated may lead to a better understanding of
the system and to suggestion for improving it, suggestions that otherwise would not be
apparent.
24 System Modeling and Simulation
4. Simulation can be used as a pedagogical device for teaching both students and practitioners
basic skills in theoretical analysis, statistical analysis, and decision-making.
5. Operational gaming has been found to be an excellent means of simulating interest and
understanding on the part of the participants, and is particularly useful in the orientation
of persons who are experienced in the subject of the game.
6. Simulations of complex systems can yield valuable insight into which variables are more
important than others in the system and how these variables interact.
7. Simulation can be used to experiment with new situations about which we have little or
no information so as to prepare for what may happen.
8. Simulation can serve as a “pre service test” to try out new policies and decision rules for
operating a system, before running the risk of experimenting of the real system.
9. When new components are introduced into a system, simulation can be used to help foresee
bottlenecks and other problems that may arise in the operation of the system.
Monte Carlo method of simulation is one of the most powerful techniques of simulation and is
discussed below.
EXERCISE
Here a,b, c, and d are given numbers, convert this model to distributed lagged model (4).
This page
intentionally left
blank
123456789
123456789
123456789
123456789
123456789
2
123456789
123456789
123456789
PROBABILITY AS USED IN
SIMULATION
Events, whose outcome can not be predicted with certainty, are called
probabilistic events. Knowledge of probability theory is a prerequisite
for the simulation of probabilistic events. Unlike the scientific experi-
ments in engineering, where outcome of results is known, in case of
random events it can not be exactly predicted. Outcome of such events
generally follow a special pattern which can be expressed in math-
ematical form called probability distribution. By probability, we mean
chances of occurrence of an event. In the present chapter, basic
concepts of probability theory are discussed in details for brushing Johann Carl Friedrich Gauss
up the knowledge of the students. In nature, number of events are (30 April 1777 – 23 February 1855)
known to be random and results for such events can not be predicted was a German mathematician and
with certainty. Let us take an example of tossing of a coin. No one scientist who contributed signifi-
cantly to many fields, including
can tell whether head or tail will be the outcome. But if we toss coin number theory, statistics, analysis,
for a large number of time, half of the time, head will be outcome differential geometry, geodesy,
and rest half of the time, tail will be the outcome. In the language electrostatics, astronomy, and
optics. He made his first ground-
of the probability, we say that the probability of coming head on top breaking mathematical discoveries
is 0.5. Similarly failure of a machine, life time of an electric bulb or while still a teenager. He completed
Disquisitiones Arithmeticae , his
impact point of a bullet on the target also follows some probability magnum opus, in 1798 at the age
distribution. Estimation of weapon delivery error is one of the im- of 21, though it would not be
portant parts of weapon performance studies. When a weapon is published until 1801. This work was
fundamental in consolidating
launched on a target from a distance, there are chances that it may number theory as a discipline and
or may not hit at a desired point. This is due to the inherent error has shaped the field to the present
in weapon release mechanism. A simple example will help to under- day.
stand weapon delivery error. If one tries to drop a stone from a height
on a small circle drawn on the ground, it is quite possible that stone may or may not fall inside the
circle. Almost same is the problem when a weapon is dropped on a target from a distance. Thus
falling of a weapon on some other point than its aim point is called weapon delivery error. Similarly
failure of a machine or one of its parts can not be predicted. In this chapter, we will study that all
the phenomena described above can be predicted, but with some probability. The problem in this
28 System Modeling and Simulation
connection is to evaluate the characteristics of a system in probabilistic or statistical terms and evaluate
its effectiveness. In section 2.4, first basic concept about various distribution functions will be given,
and then some case studies will be taken up for illustration purpose. First of all knowledge of some
basic parameters is essential so as to understand the various probability density functions. These
parameters will be discussed in sections from 2.1 to 2.3.
A B A B A
C
A
A∪ B A∩ B
Fig. 2.1: Venn diagram for union, intersection and complement of three sets.
30 System Modeling and Simulation
P(A|B) = P(A)
and
P(B|A) = P(B).
Symbol P (A|B) here means, probability of occurrence of A if B has already occurred. In other
words, if A and B are independent, then the conditional probability of A, given B is simply the individual
probability of A alone; likewise, the probability of B given A is simply the probability of B alone. This
result is after Swine and is called SWINE’S THEOREM.
We further elaborate the concept to understand this. If X is a real-valued random variable and
a is a number, then the event that X ≤ a is an event, so it makes sense to speak of its being, or not
being, independent of another event.
Two random variables X and Y are independent if and only if, for any numbers a and b the events
[X ≤ a] (the event of X being less than or equal to a) and [Y ≤ b] are independent events as defined
above. Similarly an arbitrary collection of random variables—possible more than just two of them—
is independent precisely if for any finite collection X1, ..., Xn and any finite set of numbers
a1, ..., an, the events [X1 ≤ a1], ..., [Xn ≤ an] are independent events as defined above.
2.1.7 Mutual Exclusivity
In section 2.1.6, we define intersection of event A and B as
In other words, the probability of A happening, given that B happens, is nil since A and B
cannot both happen in the same situation; likewise, the probability of B happening, given that
A happens, is also nil.
(b) Difference Rule: If occurrence of A implies occurrence of B, then P(A) ≤ P(B), and
the difference between these probabilities is the probability that B occurs and A does not:
P(B and not A) = P(BAc) = P(B) ( 1 − P ( A) )
Proof: Since every outcome in A is an outcome in B therefore A is a subset of B. Since B can
be partitioned into A and (B but not A),
P(B) = P(A) + P(BAc)
hence the result.
Now since probability of an outcome A in one trial is P(A) then probability of not occurrence
of A is 1 – P (A). Probability of not occurrence in two trials is that it does not occur in first trial
and it also does not occur in second trial, which is same as (1 – P(A))2. Similarly if we extrapolate,
probability of not occurrence of A in n trials is (1 – P(A))n. This means A does not occur in first
trial and it also does not occur in all the n trials. Therefore probability that A occurs at least once
in n trials is 1 – (1 – P(A))n. This result is very useful and will be used in next chapters.
(i ) f ( x ) ≥ 0 U|
(ii ) ∑ f ( x ) = 1V
x
|W ...(2.4a)
Here f (x) is called the Probability Density Function (PDF) of X and determines the
distribution of the random variable X. In the second equation of (2.4a) summation is over the all possible
values of x. Meaning of expression Pr(X = xj) = pj here, is that the probability that X = xj is pj . Equation
(2.4) says that distribution of random variable X is given by function f (x) which is equal to pj for x
= pj. To understand this, let us take an example of two dice. In this example, if we throw two dice
together, the probability of getting a sum two of the two faces is 1/36 (There are total 36 ways, when
two dice can fall, and sum two can come only in one way i.e., (1 + 1 = 2). Thus we can say,
f (2) = 1/36
similarly sum four can come in three ways (1 + 3, 3 + 1, 2 + 2), thus
f (4) = 3/36
and so on. Many times the PDF is in the form of a table of probabilities associated with the values
of the corresponding random variable whereas sometimes it might be expressed in some closed form,
such as
x n− x
f ( x ) = Pr( X = x ) = p (1 − p) ; x = 0, 1,..., n and 0 < p < 1 ...(2.5)
Equation (2.5) will be discussed in details in section 2.4.1. At this stage, we will assume that
f (x) is some analytic function. In the coming sections, we will discuss number of probability density
functions which are in terms of closed form functions.
µ, is the sum of the product of all the values, a random variable takes, with its probabilities assigned
at those values. Thus expected value is defined as [23]
E(X) = ΣR(x) Pr (X = x) . x
= ΣR (x) f (x) . x ...(2.6)
where R(x) denotes the range of X (Sample Space) and Σ is a symbol of summation which means,
sum of all the values over the sample space R(x). Range or sample space R is the region in which
all the values x of stochastic variable X fall. We can call expected value as the mean of values of
X in the range R. This will be clear from the exercise 2.3 given below.
Number expected (X) Probability f(x) Number expected (X) Probability f(x)
1 6
2 7
36 36
2 5
3 8
36 36
3 4
4 9
36 36
4 3
5 10
36 36
5 2
6 11
36 36
1
12
36
2.3.2 Variance
Another important quantity closely associated with every random variable is its variance. Basically,
the variance is a measure of the relative spread in the values, the random variable takes on. In particular,
the variance of a discrete random variable is defined as
Var(X) = E {X – E(X)}2
= Σ R ( X ) Pr ( X = x )[ x – E ( x )]2
2
= S R( X ) f ( x ) x – E ( x ) ...(2.6d)
This will be clear when we work out the exercise 2.4.
In Fig. 2.3, variation of f (x) versus x has been shown. Variation of f (x) is of the type of an
isosceles triangle, vertex being at x = 7, f (x) = 5.83333. This means that x = 7 is the mean value
of the variables x for which maximum spread in f (x) is 5.83333. Relation (2.6e) shows, Var(x) is
nothing but sum of the squares of deviation of x from its mean value E(x) multiplied by the probability
density function f (x) at x. Square root of Var(x) is called the standard deviation of the variable X
and is denoted by symbol σ.
d
[ f ( X )] = 0
dx
d2
[ f ( X )] < 0
dx 2
which is the condition that f (x) is maximum.
2.4.2 Median
Median divides the observations of the variable in two equal parts. Thus for a discrete or continuous
distribution of variable x, if X denotes the median, then
p(x ≤ X) = p(x ≥ X) = 0.5
Median can easily be found from Cumulative Distribution Function because at median value of
CDF = 0.5. Relation between mean, mode and median is given as,
Mean – Mode = 3(mean – median).
The relative values of mean, mode and median depends upon the shape of the probability curve.
If the probability distribution function is symmetric about the centre and is unimodal then all the three
will coincide. If the curve is skewed to the left as in the Fig. 2.4 then mode is greater than the mean,
but if it is skewed towards right then mode is greater than the mean.
A comparison of the three measures of central tendency will reveals that, the mean involves the
weightage of data according to their magnitude. The median depends only on the order, while the
mode depends only on the magnitude.
Probability as Used in Simulation 37
F(x) = ∑ f ( z)
X ≤x
...(2.7)
F(x) represents a probability that the random variable X takes on a value less than or equal to
x. We can see that in the example of rolling pair of dice only two outcomes less than or equal to
three are f (2) and f (3), thus (Example 2.1).
38 System Modeling and Simulation
∑ p( x ) = 1
R( X )
i
1
f (x) = , for a ≤ x ≤ b ...(2.8)
b–a
The mean and standard deviation of uniform distribution is given by
b
a+b
, σ = ∫ ( x – µ) f ( x) dx
2 2
µ =
2 a
b 2
1 a + b
b – a ∫a
= x – dx
2
b
1 a + b
3
= 3(b – a ) x –
2
a
Probability as Used in Simulation 39
1 b − a a−b
3 3
= –
3(b – a ) 2 2
=
(b – a )2 ...(2.8a)
12
Cumulative distribution function of uniform distribution is
x
1 x–a
F(x) = ∫ b–a dx =
b–a
...(2.8b)
a
f(x) F(x)
1(b – a)
a/(b – a)
x a x
In Fig. 2.6 we have given variation of f (x) and F(x). It can be seen that f (x) is a horizontal
straight line whereas F(x) is an inclined line. As x is increased by a, F(x) is incremented by
a/(b–a).
12, 13 and 23
This can be written as
C23 = 3
Thus combination of three numbers taken two at a time is three. We have assumed in the definition
of Cxn that combinations 23 and 32 are same. Thus Cxn denotes the number of ways k success can
occur in n independent experiments or trials and is given by
∑C
n
x p x (1– p ) n – x = 1 ...(2.12)
x=0
f ( x) = Cxn p x (1 − p) n − x ...(2.13)
Expected value and variation of binomial distribution function is given by [24,72] (Appendix-2.1)
E (X) = np
Var(X) = np(1 – p) ...(2.14)
Let us understand binomial distribution function by following example.
Example 2.5: In an examination, total thirty students appeared. If probability of passing the
examination of one student is 0.6 ( i.e., p = 0.6) , then what is the probability that none of the student
will pass and twenty students will pass.
Solution:
This problem can easily be solved by Binomial distribution. Probability of failing all the students
will be f (0) and passing of 20 students will be f (20). Thus using equation (2.13), one gets
Probability as Used in Simulation 41
This problem has a direct physical relevance. Since probability of failing one student is 0.4, this
means that probability of failing all the twenty students is remote. Thus f (0) is a very small number.
At the most 12 students are expected to fail and 18 expected to pass (E(x) = 18) out of 30 students.
In Fig. 2.6a, we have drawn probability f (x) vs. x, where x is the number of passing students. This
means in another words that probability of passing 18 students is maximum (f (18) = 0.147375) in
this case. Computer program written in C++ is given (Program 2.1) for the interest of readers. As
has been mentioned in the beginning, knowledge of numerical computation and C++ programming
is a pre-requisite for understanding the contents of this book. C++ is the most versatile and latest
language in programming. We will explain all the steps in C programming whenever they are given.
However student is advised to refer any book on C++ programming[3].
Program 2.1: Binomial Distribution Program
#include <iostream.h>//Header files to be included. Note that comment command
can be used anywhere.
#include <math.h>
double n;
42 System Modeling and Simulation
cout<<“c_n_x = “<<c_n_x<<endl;
// Writes on console the value of C nx . Endl means end of line .
}
void comp::fx()
{
/*When ever we write a function defined in a class we write class name::
function name. Symbol :: is called scope resolution operator. Function fx()
computes Binomial probability distribution function
f ( x) = Cxn p x (1– p) n – x. .*/
double fx;
fx= c_n_x* (pow ( p ,x )*( pow((1– p),( n– x))));
// Function pow(p, x) means p x.
cout<<“fx =” <<fx <<endl;
}
double comp::fact(double x)
{
/*This function computes the factorial of x. First we define fact=1 as initial
value. Function for (int j=x ; j > 0; j--) says multiply fact with j and store to
fact in each loop, where j varies from x to 1. j–means decrease the value of j by one
every time control enters the for () loop*/
double fact = 1.0 ;
for (int j = x ; j >0 ; j --)
{
fact *= j ;
}
return fact;//Returns the final value of factorial x.
}
void main()
{
/* Object obj of class comp has been declared.
44 System Modeling and Simulation
comp obj;
obj.pre_input();/* Will call function pre_input()defined in class
comp*/
for(int i=0;i<=n;++i)
{
//for() loop is executed n times.
obj.input();
obj.C_n_x();
obj.fx();
}}
e– λ λx
= , where x = 0, 1, 2,... ...(2.15)
x!
and λ > 0.
Here λ is the expected value of the Poisson’s distribution. Poisson’s distribution is used in the case
where out of large number of trials, probability of success of an event is quite low. For example, if a
shell explodes in the near vicinity of an aircraft, probability of hitting only two or three fragments to a
vital part, out of large number of fragments is too small and is calculated by using Poisson’s distribution
[54]. It can be proved that sum of all the probabilities when x varies from zero to infinity is unity i.e.,
n
e– λ λ x
∑
x=0 x!
= 1 when n → ∞
Proof:
n
e– λ λ x n
λx
∑ x!
= e ∑
–λ
= e–λ. eλ = 1
x=0 x = 0 x!
n
λx
where ∑
x = 0 x!
= 1
1. Sime’on Denis Poisson reported this distribution for the first time in his book “Recherches sur la probabilite des
judgements (1937)” (Researches on the probability of opinions) [20].
Probability as Used in Simulation 45
It can be shown that, if we put p = λ/n where n is very large, binomial distribution reduces
to Poisson distribution i.e., as n tends to infinity.
Expression for Binomial distribution function is
n( n − 1).( n − 2)...( n − x + 1) x
= λ (1 – λ/n)n– x
n x . x!
1 2 x − 1
1 − . 1 − ... 1 −
n n n λx ( 1 – λ /n ) n − x
=
x!
If we now let n → ∞, we get
1 . 2 ... x − 1
1 − 1 − 1 − →1
n n n
and
n–x λ
λ λ n / λ λ
−x
1– = 1 − 1 − → e−λ
n n n
λ n / λ λ
−x
e −λ λ x x = 0,1, 2,...;
= ; ...(2.16)
x! λ>0
Expected value and variance for Poisson random variable are
E(X) = λ
Var (X) = λ ...(2.17)
Below, we will demonstrate by a numerical example the validity of equation (2.16).
In the Table 2.3, we have given a comparison of both the distributions for n = 100 and λ = 1.0
We can easily see from the table that under the given conditions, both the distributions give similar
results. Thus it is proved that Poisson’s distribution is a special case of binomial distribution.
x 0 1 2 3 5
Binomial .366 .3697 .1849 .0610 .0029
Draw the probability distribution histogram for X = the age of a randomly chosen college student.
Solution:
Summing the entries in the bottom row, we see that the total number of students in 1980
was 12.356 million. We can therefore convert all the data in the table to probabilities dividing
by this total.
0.39
0.22
0.16
0.13
0.10
We can define continuous random variable in another way too. A random variable X and the
corresponding distribution of random variable is to be continuous if the corresponding cumulative
distribution function F (x) = P(X ≤ x) can be represented by an integral form.
x
F (x) = ∫ f ( x) dx ...(2.18)
−∞
where the integrand f(x) is continuous everywhere except some of the finite values of x. The integrand
f(x) is called the probability density function or density of the distribution. In equation (2.18) we see
that F'(x) = f (x) for every x at which f (x) is continuous, where F′ (x) is the derivative of F(x) with
respect to x. In this sense the density is the derivative of the cumulative distribution function F(x).
Since sum of all the probabilities in the sample space is equal to unity, from equation (2.18) one gets,
+∞
∫
–∞
f ( x)dx = 1 ...(2.19)
Thus probability P(a < X ≤ b) is the area under the curve with function f (x) and between the
lines x = a and x = b (see Fig. 2.7).
25
20
15
f(x)
10
x=a x=b
0
Fig. 2.7: Probability function P(a < x ≤ b), is area under the curve f(x) and lines x = a and x = b.
48 System Modeling and Simulation
Obviously for any fixed a and b (> a) the probabilities corresponding to the intervals a < X ≤ b,
a ≤ X < b, a ≤ X ≤ b, a < X < b are all same. This is different from a situation in the case of a
discrete distribution. Since the probabilities are non-negative and (2.18) holds for every interval, we
must have f (x) > 0 for all x.
e– x / λ
,t ≥ 0
f(x) = λ ...(2.21)
0, elsewhere
The graphical representation of density function is shown in Fig. 2.8a. The mean and variance
of exponential distribution has given by
2
E(X) = µ = λ and Var(X) = λ ...(2.22)
If variable x is time t, then exponential density function sometimes is written as,
1 –t / τ
e , t ≥0
f (x) = τ ...(2.23a)
0, elsewhere
where τ = λ is the mean of the distribution, also called mean time of arrival.
Cumulative probability distribution function of exponential density distribution function is given
t
1 –t / τ
τ ∫0
F(t) = e dt = 1– e – t / τ ...(2.23b)
1 1
f (x)
F(x)
0
t/ t/
(a) Exponential density function (b) Exponential distribution function
The exponential distribution has been used to model inter arrival times when arrivals are completely
random and to model service times which are highly variable in the theory of queuing. Queuing theory
will be discussed in chapter seven. This distribution is also used to model the life time of a component
that fails catastrophically, such as light bulb. In this case λ is the failure rate.
x α −1 e − x / β x ≥ 0
f (x) = ...(2.24)
Γ (α )β α α , β ≥ 0
where α, β are positive parameters and
∞
∫x
α −1 − x
Γ (α ) = e dx for all a > 0 ...(2.24a)
0
This result is independent of λ (since x = λ). Thus we observe that probability that bulb will
last more than its mean value is 0.368, for any value of λ.
Probability that electric bulb will last between 2000 and 3000 hours is,
P (2 ≤ X ≤ 3) = F (3) − F (2)
−3/ 3 −2 / 3
= (1 − e ) − (1 − e )
= −0.368 + 0.513 = 0.145
Let us now discuss one of the most important property of exponential density distribution i.e.,
it is “memory less”. That means, the probability that a component lives at least for t + s hours, given
that it has already lived for s hours, is same as the component at least lives for t hours, if life of
component is exponentially distributed. Which means that for all s ≥ 0 and t ≥ 0,
P ( X > s + t / X > s ) = P( X > t ) ...(2.25)
In equation (2.25), X represents the life of a component (a bulb say) and assume that X is
exponentially distributed. Left hand side of the equation (2.25) states that the probability that the
component lives for at least s + t hours, given that it has already lived for s hours, and right hand
side states that probability that it lives for at least t hours. If the component is alive at time s
(if X > s) then the distribution of remaining amount of time that it survives, namely X – s, is the
same as the original distribution of a new component. That is, the component does not “remember”
that it has already been in use for a time s. A used component is as good as new.
That equation (2.25) holds is shown by examining the conditional probability
P ( X > t + s)
P( X > s + t / X > s) = ...(2.26)
P( X > s)
From equation (2.22b), we know that F(t) is the probability that X is less than or equal to t.
Thus probability that component lives at least for t hours is equal to 1 – F(t).
Substituting 1 – F(t) = e–t/λ in the right hand side of (2.26), we get
e− ( s +t ) / λ
P( X > s + t / X > s) = = e−t / λ ...(2.26a)
e− s / λ
which is equal to P(X < t).
Example 2.8: Find the probability that industrial lamp in example 2.7 will last for another 1000
hours, given that it is operating after 2500 hours.
Solution: We use equations (2.22) and (2.22a) and get,
P( X > 3.5 / X > 2.5) = P( X > 1) = e−1/ 3 = 0.717
This logically is not correct as time increases, probability of its functioning goes on increasing.
For example probability of surviving for 6000 hours is nothing but e–2.0 = 0.135134.
2.7.2 Erlang Density Function
It has been mentioned while discussing Gamma density function that when α is an integer it becomes
Erlang density function. A. K. Erlang was a Danish engineer responsible for the development of queuing
λ
theory. Let α = k, and β = in equation (2.24), where k is an integer. Equation (2.24) becomes
k
x k −1e − kx / λ x ≥ 0
f (x) = (k / λ ) k ...(2.27)
(k − 1)! k , λ ≥ 0
Probability as Used in Simulation 51
1 – x/λ
f (x) = e , x ≥ 0, λ ≥ 0
λ
which is nothing but exponential distribution.
z
+•
In (2.27a) f (xj) is the probability function of random variable X. The mean is also known as
mathematical expectation of X and is denoted by E(X). In case of continuous random parameter mean
is nothing but integral over all the values of X.
A distribution is said to be symmetric with respect to a number c if for every real x, if
f(c + x) = f(c – x).
The variance of distribution is denoted by σ2 and is defined by second moments as follows:
(a) σ2 = ∑ j (xj – µ)2 f(xj) (discrete distribution) …(2.28)
z
+•
Solution: In the table xi , yi are experimental values and y are those obtained from the fitted
expression. Variance of y thus is
52 System Modeling and Simulation
( y − yi )2
s = ∑ n
thus
.0828
∴ σ = s = = 0.0909
10
Thus standard deviation of experimental data from a fitted curve is 9.09%.
1
f (x) = exp{– ( x – µ) 2 / σ 2 }, σ > 0 ...(2.29)
σ 2π
is called the normal distribution or Gaussian distribution. A random variable having this distribution
is said to be normal or normally distributed variable. This function is like an inverted bell shape being
symmetric about point x = µ and is shown in Fig. 2.9.
This distribution is very important, because many random variable of practical interest are normal
or approximately normal or can be transformed into normal random variables in a relatively simple
fashion. In equation (2.29) µ is the mean and σ is the standard deviation of the distribution.
Probability as Used in Simulation 53
Φ(z) = e–u /2
du ...(2.32)
2π –∞
which is the distribution function with mean equal to zero and variance equal to one. This integral
x−µ du 1
is obtained from integral (2.31) by substituting, = u, = , and limits for integration varying
σ dx σ
x−µ
from −∞ to z= .
σ
Therefore equation (2.32) becomes
( x − µ)
1
∫
2
F ( x) = σ
e− u / 2 du ...(2.33)
2π −∞
x−µ
z =
σ
x−µ
F ( x) = Φ
σ
Therefore combining eqs. (2.31) and (2.33) one gets
b − µ a − µ
P( a < X ≤ b) = F (b) − F ( a) = Φ − Φ ...(2.34)
σ σ
in particular, when
a = µ−σ
and b = µ+σ
we get P (µ − σ < z ≤ µ + σ) = Φ (1) – Φ (–1)
54 System Modeling and Simulation
Integral in equation (2.32) have been integrated numerically by using Simpson’s method [4] and
results are shown as follows
(a) P (µ − σ < x ≤ µ + σ ) ≈ 68%
(b) P (µ − 2σ < x ≤ µ + 2σ ) ≈ 95.5%
(c) P (µ − 3σ < x ≤ µ + 3σ) ≈ 99.7%
A computer program for normal distribution function written in C++ language by Simpson’s
method is given below.
sum += 2.0*func(x);
}
area = h*sum/3.0;
diff = fabs(prev_result - area);
prev_result = area;
cout<<n<<’\t’<<h<<’\t’<<area<<’\t’<<diff<<’\t’<<acc;
/*Here ‘/t’ means give a space between the value of parameters n,
h,area,diff,acc.*/
n +=2;
}while (diff > acc);
if (diff <= acc )
{
cout<<“\n \n”);
cout<<“ integral of the function e(X^2) from a to b is \n”;
cout<<area;
getch();
}
return;
}
In the following paragraph we demonstrate an experiment on computer which demonstrates the
normal distribution.
In the simulation example numbers obtained in different circles are slightly different from
theoretical values. This is because the total number of trials are not large (n = 2000). If n is increased,
closer results will be obtained.
1 n
sx2 = ∑ ( xi − x )2
n i =1
...(2.35)
1 n
Sx = ∑ ( xi – x )2
n i =1
The standard deviation in this case is not the true standard deviation. It always has some error
and is different from true standard deviation, which will be denoted by σx. We call this as biased
standard deviation. But if the sample size is very large then this error reduces to zero. Thus true standard
deviation is defined as,
1 n
sx = Lim
n →∞ n
∑
i =1
( xi – x )2
s = s x2 + sy2 = 2.898
If we consider a very large number of data points, then standard deviations Sx and Sy would
approach to their true values σx and σx respectively.
It is known that sample variance is generally computed as [72]
1 n
Sx = ∑ ( xi – x )2
n − 1 i =1
...(2.36)
which is based on (n– 1) degrees of freedom—one degree of freedom being used in the calculation
of the sample mean x as an estimation of the population mean. This is for large values of n and
is equal to σ 2x . In fact both the sample variances converge to same value when n is large
(see Appendix 4.2 : Sampling distribution of means).
The sample means and standard deviation describe only the location of the centre of data points
and their dispersion. These two parameters do not however describe the characteristics of the overall
distributions of data points. In case data points are hits on a target due to a weapon, we have to make
58 System Modeling and Simulation
some assumptions of reasonable practical value, in estimating probabilities of hitting. It is widely assumed
and also verified sufficiently that the distribution of rounds is approximately normal or Gaussian in
character. Thus, the probability density function of rounds in x-direction may be described by
2
1 (x − x )
f ( x) = exp − 2 ...(2.37)
2πσ2x 2σ x
and that of the dispersion in y-direction by
2
1 ( y – y )
f ( y) = exp 2 ...(2.38)
2π σ 2y 2σ y
1 x 2 y 2
f ( x, y ) = exp − 2 − 2 ...(2.38a)
2πσ x σ y 2σ x 2σ y
If σx ≠ σy, we call this as non-circular and if σx= σy then it is circular bi-variate normal distribution
and distribution is given by,
1 x2 + y 2
f ( x, y ) = exp − ...(2.38b)
2π σ 2 2σ 2
which shows that (x, y) is the point normally distributed around origin.
2.10 CIRCULAR PROBABLE ERROR (CEP) AND THE PROBABLE ERROR (PE)
In this section we discuss various errors in case of weapon delivery which will be of interest to defence
scientists. A measure of dispersion generally used to describe weapon delivery accuracy is the Circular
Error Probable (CEP). Parameter CEP is the basic parameter for determining the error in weapon
performance and is defined as the radius of the circle about the true mean point with respect to aim
point; which includes 50% of the hits, considering a very large number of rounds fired onto the target
area under stable firing conditions.
Given the circular normal density function (equation (2.38b))
1
f ( x, y ) = exp[ − ( x 2 + y 2 ) / 2σ 2 ] ...(2.39)
2πσ 2
we integrate this function over a circular target with the centre at the origin and equate the result
to one-half
where radius R0.50 is radius of circle so that 50% of the hits fall inside it.
By making the transformations of variables
x = r cosθ, y = r sinθ, 0 ≤ θ ≤ 2π
in the integral (2.40) one gets,
R0.5 2 π
1
∫ ∫ d θ exp [ − r
2
= / 2σ 2 ] r dr ...(2.41)
2πσ 2 0 0
= 1 − exp − R0.5
2
/ 2σ2 = 0.5
or
exp − R0.5
2
/ (2σ 2 ) = 0.5
1 x2 + y 2
f (x, y) = 2 exp − ...(2.42)
2π σ 2σ 2
Here σ is the standard deviation of weapon. Then the chance of a round hitting the circular target
is simply (see section 2.10, equation (2.41))
Equation (2.43) only tells whether centre of lethal area of a particular weapon falls on a target
whose radius is R or not and if it falls what is the probability. (This probability is also the chance
that chi-square with two degrees of freedom is χ 2 (2) and it does not exceed R2/σ2 ). To understand
relation (2.43), let us consider the following example.
Example 2.10: Fire from an enemy bunker is holding up the advance of friendly troops. If an
artillery with a warhead damage radius of 30m and delivery CEP of 20m is used what is the chance
of destroying the bunker?
Here it is essential to tell that standard deviation σ is function of distance of aim point from the
launch point of the weapon (Range) and is written as σ = rθ , where θ is the dispersion of weapon
in radians and r is the range. But if range is known beforehand, CEP can be given in terms of distance
in place of angle.
Solution:
CEP = 20 = 1.1774 σ
Therefore σ = 17,
The chance of killing the point target (bunker is a point target) may be found from the chance
of a round falling on or within the radius of 30m from it. Thus
302
p ( h) = 1 − exp − = 0.789774
2 × 17 × 17
It can be seen from the above equation that chance of hitting the target is 0.79 i.e., there are 79%
chances that bunker will be destroyed.
Note that in this example we have deduced that with 79% probability, bunker will be destroyed.
This logic does not seem to be correct, because target is a bunker, which is not a levelled target.
Due to the depth of the bunker, it may not be possible to destroy it. In fact equation (2.34) only
gives us that target is covered by the weapon.
Probability as Used in Simulation 61
EXERCISE
1
x , for 0 < x < 2
f (x) = 2
0, otherwise
Plot the associated probability distribution histogram, and use it to evaluate (or estimate)
the following:
(a) P(0 ≤ X ≤ 4)
(b) P(X ≥ 4)
(c) P(2 ≤ X ≤ 3.5)
(d ) P(X = 4)
11. Generate three random variates from an exponential distribution having mean value 8.
Hint: A variate of exponential distribution is given as,
1 1
y = − ln(1 − r ) or − ln( r )
µ µ
where 1/µ is the mean of the distribution and r is the uniform random number.
12. There are 15 equally reliable semiautomatic machines in a manufacturing shop. Probability
of breakdown per day is 0.15. Generate the number of break down for next seven days.
Determine the mean and variance of the generated observations. Compare with the
theoretical values of mean and variance.
13. The life time, in years, of a satellite placed in orbit is given by the following exponential
distribution function,
0.4e −0.4 x , x ≥ 0
f (x) =
0, otherwise
(a) What is the probability of life of satellite being more than five years?
(b) What is the probability of life of satellite being between 3 and 6 years?
14. The distribution function for a random variable x is:
3 − e − x , x ≥ 0
F (x) =
0 , x <0
Find:
(a) Probability density function
(b) P(x > z)
(c) Probability P(–3 < x ≤ 4). (PTU, 2002)
15. Give expressions for Binomial, Poisson and Normal distributions. Under what conditions
Binomial distribution is approximated by Poisson distribution. (PTU, 2002)
16. 5000 students participated in a certain test yielding a result that follows the normal
distribution with mean of 65 points and standard deviation of 10 points.
(a) Find the probability of a certain student marking more than 75 points and less than
85 points inclusive.
(b) A student needs more than what point to be positioned within top 5% of the participants
in this test?
(c) A student with more than what point can be positioned within top 100 students?
Probability as Used in Simulation 63
APPENDIX 2.1
∑C
x=0
n
x p x q n− x = 1
Then
n
E (X) = ∑∑ x f ( x, p, n)
x=0
n
n!
= ∑x x=0 x ! (n − x)!
p x qn− x
n
(n − 1)!
= np ∑
x =1 ( x − 1)! ( n − x)!
p x −1q n − x
since the value of term with x = 0 is zero. Let s = x – 1 in the above sum. Thus
n −1
( n − 1)!
E (X) = np ∑
s=0 s ! (n − 1 − s )!
p s q n −1− s
Now since
n −1
(n − 1)!
∑s =0 s !(n − 1 − s )!
p s q n −1− s = ( p + q)n−1 = 1
we get
E(X) = np
x=0
64 System Modeling and Simulation
n
n!
= ∑ x . x!(n − x)! p q
x=0
2 x n− x
n
x ( n − 1)!
= np ∑ ( x −1)! ( n − x)! p
x =1
x −1
q n− x
s=0
n −1
= np ∑ (s + 1)
s=0
f ( s, n − 1, p )
n –1 n –1
n –1
= ∑ f (s, n − 1, p)
x=0
= (n − 1) p + 1
= np + q
Thus E ( X 2 ) = np (np + q )
and we get
2 2 2 2
Var( X ) = E ( X )2 − ( E ( X ))2 = n p + npq − n p = npq
12345678
12345678
12345678
12345678
12345678
3
12345678
12345678
12345678
AN AIRCRAFT SURVIVABILITY
ANALYSIS
Aircraft survivability analysis is an important study in the aircraft industry. It means, an aircraft when
under enemy attack is capable of surviving or not. This study is useful for estimating the damage caused
to the aircraft, as well as deciding the basic parameters of the aircraft, while under design stage. In
chapter two we have studied probability distribution of discrete random events. Also various laws of
probability have also been studied. In this chapter we will apply these laws for making a model of aircraft
survivability. In chapter two we explained how circular normal distribution can be applied to calculate
single shot hit probability. In this model these concepts will be further elaborated. Due to the complex
nature of the aerial targets, their damage analysis is an independent field and it is necessary to deal it in
a different chapter. Under aerial targets, various airborne vehicles such as aircraft, helicopter, Remotely
Piloted Vehicles (RPVs) etc. can be categorised. There is an important term called ‘denial criteria’ in
damage assessment problems. Denial criteria means how a target will be damaged, partially, completely
and will be unserviceable. It has been observed from war data that even if an aircraft is hit and had
several bullet holes, has come back from the battle field. Therefore the criterion that aircraft or any
other target under attack is completely destroyed is very important and requires a separate study. A denial
criterion for the aircraft is quite different from the ground targets and requires a detailed study. Aircraft
consists of hundreds of vital parts and kill (totally damaged) of either of them can lead to aircraft’s
malfunctioning, which can lead to its kill after sometime of operation. It is quite important to study each
and every vital part of the airborne vehicle critically, the energy required, and type of mechanism needed
for its kill, all of which matter in this study. In this chapter we will discuss the basic mathematics required
for modelling the aircraft vulnerability. Since aircraft is a moving object, and it has various capabilities to
avoid enemy attack, the study of its survivability against air defence system is called, aircraft combat
survivability study. Aircraft combat survivability is defined as the capability of an aircraft to avoid
and/or withstand a man made hostile environment. This study is of quite importance during the design
stage of aircraft, as well as under its operational stage. In the present chapter a general concept of
susceptibility of aircraft, its vital parts and their kill mechanism alongwith various kill criteria will be
discussed. A dynamic model of aircraft vulnerability will be taken up in chapter six. In this chapter only
66 System Modeling and Simulation
a static model of aircraft survivability, using probabilistic concepts learnt so far will be used to build the
model. Here it will be assumed that projection of aircraft as well as its vital parts on a plane has been
provided as inputs. In chapter six, we will study with the help of mathematical models, how to obtain
projection of aircraft on a given plane. An aircraft has hundreds of vital parts but for the sake of simplicity,
only three vital parts are considered in this model. More elaborate model will be discussed in chapter six.
1. Blast and fragments strike the aircraft Yes How many fragments hit the aircraft and
where do they hit?
2. Missile warhead detonates within the Yes Can the onboard Electronic Counter Measure
range? (ECM) suite inhibit the functioning of the
proximity fuse?
Contd...
An Aircraft Survivability Analysis 67
3. Radar proximity fuse detects aircraft Yes Will chaff decoy the fuse?
4. Missile propelled and guided to vicinity Yes Can the target aircraft outmanoeuvre the
of aircraft. missile?
5. Missile guidance system functions in Yes Are i.r. (infra-red) flares effective decoys?
flight.
6. Missile guidance system locked on to Yes Is the engine’s infra-red suppresser effective
engine (infra-red) i.r radiation. in preventing lock on?
7. Target’s engines within missiles field of Yes Are the engine hot parts shielded?
view?
8. Enemy fighter manoeuvres to put target Yes Does the enemy fighter have a performance
into field of view and within maximum edge? Does the target A/C have an offen-
range. sive capability against the enemy fighter?
9. Target aircraft designated to enemy Yes Does the on-board or stand off ECM suite
fighter and fighter launched. have a communications and jamming
capability?
10. Fighter available to launch against the Yes Are there any supporting forces to destroy
target. the enemy fighter on the ground?
11. Enemy ‘C3’ net function properly. Yes Does the stand-off ECM suite have a
communications jamming capability?
(b) Tactics: Tactics is another parameter which depends on pilot’s skill and helps in reducing
the susceptibility. To avoid detection by the enemy air defence system, pilot can hide his
aircraft behind the terrains. This technique is called terrain masking. Also an aircraft sortie
consisting of escorts aircraft to suppress enemy air defence, also helps in reducing the
susceptibility. It can carry various survivability equipment, such as electronic jamming
devices to jam the enemy detection system or weapon to destroy it. Designing aircraft in
such a way that it has minimum Radar Cross Section (RCS) is an art of the day. Radar
cross section of a target will be discussed in section 3.3. Stealth aircraft developed in USA
have minimum radar cross section and are difficult to be detected. These points are
explained by the following illustration of an aircraft on a mission to drop paratroopers in
enemy’s territory.
Example: Consider a case of a transport aircraft attempting to deliver troops on a bright sunny
day to a location near the FEBA (Forward Edge of Battle Area).
– It drops down into a valley to take advantage of terrain masking,
– A self propelled RADAR - directed Anti Aircraft Artillery (AAA) system detects the aircraft
with the scanning radar through its optical tracker.
– Meanwhile RWR (Radar Warning Receiver) in aircraft detects scanning signals from AAA
(Anti Aircraft Artillery) radar and alerts pilot for the type, location and status of the THREAT.
– Pilot ejects chaff, attempt to break the lock of tracking radar by manoeuvring and looks
for hide out (terrain or vegetation).
– Observer on the AAA platform, watching the aircraft, redirects radar towards it.
– After short time, when a fire control solution is obtained ground system fires at aircraft.
68 System Modeling and Simulation
The MODELING and quantification of individual events and elements in such an encounter is
referred as a SUSCEPTIBILITY ASSESSMENT. Obviously, there are many diverse factors that
influence susceptibility, many of which are difficult to model and to quantify. In order to determine
which of the factors or events and elements are the most important and which ones are of lesser
importance, an Essential Elements Analysis (EEA) should be conducted. To illustrate the EEA, consider
an encounter between a friendly strike aircraft and an enemy interceptor carrying an air-to air, infrared
homing missile. Details of the analysis is given in Table 3.1.
Let us assume that a radar transmits a pulse with power Pt , towards a target which is at a distance
R from it. Assuming that the radar transmitting antenna transmits power isotropically in all the directions,
with spherical symmetry (area of the surface of enveloping sphere being 4πR2), then power flux per
unit area at range R will be
Pt
Power density =
4πR 2
If the radar antenna has a gain G , then the power density at the target will be multiplied by G.
Now assume that target reflects back all the power intercepted by its effective area, and the reflected
pattern is also isotropic. If the reflected area of the isotropic target is σ, which is called the
radar cross-section, then the power reflected isotropically from the target is given by,
PGt
Reflected echo signal = σ
4 πR 2
Radar cross-section σ has units of area and can be measured experimentally. Since the power
is reflected isotropically, the reflected power density received by the antenna is
Also it is known that relation between gain G, and effective area A is [32]
Gλ 2
A = ...(3.2)
4π
Substituting A in (3.1) from (3.2), one gets
2 2
PG
t λ σ
PR = 3 4 ...(3.3)
(4π) R
In the derivation of radar’s equation (3.3), it has been assumed that reflections from the target are
isotropic, which is only an assumptions and not happens in actual situation. But still these equations
can be used. We replace isotropic target for original target and change the area of isotropic target in
such a manner that it gives same reflection echo, as original target gives. The area of the equivalent
isotropic target calculated in this way is called the radar cross-section of the target, which off course
is different from the actual cross-section of the target. Here in radar equation, noise factor has not been
considered. When effect of signal to noise ratio is also taken into account, equation (3.3) becomes [32,71]
1/ 4
PG t
2
λσF 4
Rmax = ...(3.3a)
(4π) Ls ξ m / n NLa
3
where
Ls = radar system losses (Ls ≥ 1),
F = relative electrical field strength of echo at the receiving antenna,
N = noise of internal system or detection process,
La = signal and echo loss due to radar transmission in the atmosphere (≥ 1).
70 System Modeling and Simulation
ξmin = (S/N) = signal to noise ratio associated with a specified probability of detection.
For the derivation of this equation, reader is advised to refer some book on radar [71].
For the s-th scan, there will be some probability of detection Pd (s) based upon actual signal to
noise ratio, which can be determined from the radar equations and cross-section, both of which can
vary with each scan. Thus Pd (s) given by equation (3.3b) is generally a function of time. If we denote
this by Pdt then probability that target has been detected after s scans is given by
s
Pd ( s) = 1 − ∏ (1 − Pdt )
t =1
B - kill: Damage caused to an aircraft that it falls out of manned control in 30 minutes,
after being hit.
These kills will be discussed in greater details in sixth chapter.
2. Aircraft description: Another factor is the assembly of the technical and functional description
of the vital parts of the aircraft. That is, to list all the vital parts, their location, size, material
function and criteria of their being killed. The functional description should define the functions
provided by each component including redundancies. Vital part, sometime called critical
component is defined as a component which, if either damaged or killed, would yield a defined
or definable kill level of aircraft. For example, pilot is a critical component for KK-type of
kill of an aircraft. If pilot is killed by a fragment or projectile, there is no hope of aircraft’s
survivability. Similarly engine in a single engine aircraft is a critical component for the A type
kill of an aircraft. But for twin engine, it is not that much critical component. It is possible
that both the engines are damaged by fragment hits. In that case even redundant engine become
critical components. To find the kill of aircraft due to the damage of it’s critical component
a Critical Component Analysis (CCA), in tabular form is to be conducted. If a particular critical
part is killed, what will be the effect on overall performance of the aircraft. This will involve
listing of all the critical components and their roll in the performance of the aircraft. This
process is called critical component analysis. Thus for critical component analysis, the first
step is to identify the flight and mission essential functions that the aircraft must perform
in order to continue to fly and to accomplish mission.
where Api , Pk/hi are respectively the projected area of i-th component on a given plane and probability
of kill of i-th component, subject to a hit on it. In the present model Pk/hi will be assumed to be
given input. Here after, all the parameters ending with capital subscripts denote values for aircraft
and that ending with lower characters give values for component. Kill probability of i-th vital component
given a random hit on the aircraft is given as,
Pk/Hi = Ph/Hi Pk/hi ...(3.6)
where
Pk/Hi = probability of kill of i-th vital part, subject to a hit on aircraft,
Ph/Hi = probability of hit on the i-th vital part, subject to a hit on the aircraft,
Pk/hi = probability of kill of i-th vital part, subject to a hit on it.
Equation(3.6) can be read as, probability of kill of i-th part of aircraft subject to a hit on aircraft
is equal to probability of hit on the i-th part of aircraft subject to a hit on aircraft and probability
of its kill subject to hit on it. Here projectile is aimed at the geometric centre of the aircraft. Probability
of hit of i-th vital part is given for simplicity, by the ratio of the projected area of the vital part and
aircraft respectively i.e.,
Ph/Hi = Api /AP ...(3.7)
where Api , AP are respectively are projected areas of vital part and aircraft, on the given plane. It
is important to mention here that in order to find projected area of a component or of full aircraft,
large number of transformations is needed, which have not been explained here. This algorithm works,
when projected areas on a plane have been provided as part of data. In chapter six, a detailed algorithm
for determining the projected areas on any plane is given. Since in actual dynamic conditions, aircraft,
projectile and its fragments after explosion, all are moving with respect to each other, a series of
projections are to be evaluated for determining a hit.
From equation (3.6) with the help of eqs. (3.5) and (3.7), one gets the probability of kill of
i-th vital part, subject to a random hit on the aircraft as,
Api Avi Avi
Pk/Hi = = ...(3.8)
AP Api AP
which is the ratio of vulnerable area of i-th vital part to projected area of whole aircraft.
It is to be noted in equation (3.5) that Avi = Api Pk/hi , which assumes that Pk/hi is the given input.
From equations (3.5) and (3.8) one gets
Api Pk / Hi
=
AP Pk / h i
In the sixth chapter, we will show that this assumption is too vague, as given a hit on the vital
part does not necessarily lead to a kill. Kill depends on the type of weapon and energy released by
it on the aircraft surface. Some of the vital parts are hidden behind other non-vital parts. A projectile
has to penetrate all these part and reach the vital part with sufficient critical energy required to kill
74 System Modeling and Simulation
it. Now we define an other term, that is probability of survival of aircraft. Probability of survival of
aircraft PS/H is defined as one minus probability of its kill, thus
PS/H = 1 – PK/H ...(3.9)
Let us assume the aircraft has n vital parts. It is known that if one of the parts is killed, aircraft cannot
fly. Thus for an aircraft to survive an enemy threat, all its critical components have to survive i.e.,
n
where
PS/H = probability of survival of aircraft subject to a random hit on it.
PS/Hi = probability of survival of i-th vital component subject to a hit on the aircraft.
For the sake of simplicity, we assume that aircraft has only three critical (vital) components i.e.,
engine, pilot and fuel tank. Equation (3.10) with the help eq. (3.9) in this case reduces to,
Throughout this chapter, we will consider an aircraft with three vital components. Below we will
illustrate these equations with numerical examples.
Case-1: To enhance the survivability of aircraft, generally more vital components are concealed
behind less vital components to protect them. In this section we will assume that vital components
are open to attack and then conceal them behind each other to see the effect of overlapping. During
flight of aircraft vis a vis projectile, positions of various components keep on changing. Let us first
consider a case when vital components of the aircraft are non-overlapping and non-redundant. By
non-overlapping we mean, while projecting aircraft on a plane, projection of two components of aircraft
are not overlapping. This condition is only a hypothetical condition as, if in one orientation two
components are not overlapping, does not mean that they will not overlap in other orientation. Actually
plane on which projection is taken is a plane, normal to the flight path of the projectile. This plane
is called normal plane (N-plane). The condition of overlapping of parts varies from fragment to
fragment. Now if components are non-overlapping, then more than one component cannot be killed
by a single penetrator i.e., only one Pk/Hi will be non-zero, other two will be zero. Here Pk/Hi = 0
means, i-th part is not killed, subject to hit on the aircraft.
Thus all the products of kill probability are zero, thus equation (3.12) reduce to,
3
PS/H = 1 – ∑ Pk / H i
i =1
let us try to understand equation (3.12) by an example. This situation is explained in Table 3.2, where
projected areas of aircraft and their vital components have been provided in meters, as a given input.
In Table 3.2, we have used the relations (3.7) and (3.8), i.e.,
Avi
Avi = Api . Pk/hi ; = Pk/Hi ...(3.12a)
Ap
In this table, in third column Pk/hi, have been assumed based on practical experience. For
example, if pilot is killed, aircraft cannot fly and thus kill of pilot is taken as 100% kill of aircraft,
that is Pk/hi = 1.0. Same way other probabilities have been provided. AP in the last row is equal to
the total projected area of aircraft which is again a given input. First column of the table gives
projected areas of three vital parts on a given plane, and column two gives their kill probabilities
subject to a hit. Here it is to be noted that a hit does not always result in total kill. There are various
other factors, such as total kinetic energy of the projectile. Column three and four are obtained
by equation (3.12a). Total kill probability of aircraft in this case is sum of the kill of three vital
parts and is 17.33% and total vulnerable area is 5.2 metre square.
Therefore, probability that aircraft’s loss is due to fire or/and fuel ingestion is = 1 – (1 – 0.3) ×
(1 – 0.1) = 1 – 0.63 = 0.37. Thus probability of fuel tank has increased from 0.3 to 0.37 due to
multiple failure modes (Table 3.3). Thus we see that in order to reduce the kill of the aircraft, engine
and fuel tank should have sufficient separation. It is to be noted that this kill probability is not the
simple sum of two probabilities. Reason for this is that due to single hit both the events may or may
not occur. Therefore the probability of not occurring of both the events is first calculated and then
probability of occurrence of at least one of the event is calculated.
Thus in Table 3.3, kill probability due to fire in fuel tank increases from 0.3 to 0.37 and total
kill probability of aircraft increases from 0.1733 to 0.1873 (Table 3.3). Thus here we see the benefit
of simple modelling, which gives direct conclusion that engine and fuel tank in an aircraft should be
at safe distance.
From Table 3.4 it can be seen that overall probability of kill reduces, due to overlapping of two
components. Due to overlapping of vital components it is quite possible that component ahead is
partially pierced and component behind it is not damaged. This will further reduce vulnerability. It
is a normal practice in modern aircraft to hide more vital components behind less vital components
to reduce its vulnerability. In this section we have assumed that there was no fire in engine due to
ingestion of fuel on the hot engine and thus engine did not received damage due to fire. In the next
section, we will extend this model by assuming that engine catches fire due to ingestion of fuel.
Thus vulnerability is enhanced due to secondary kill mode and it further increases due to fire.
Results of this section suggest that engine and fuel tank should always be kept apart. This is due
to the fact that engine becomes very hot and is always prone to fire if fuel ingestion is there. Now
a days aircraft with self sealing materials are available, which can reduce this type of risk.
3.6.2 Redundant Components with no Overlap
So far we have assumed that aircraft has only non-redundant components. It is quite possible that
aircraft has redundant components. For illustration purpose, consider an aircraft with two engines.
This means if one engine fails, aircraft still can fly with other engine. The kill expression for the aircraft
model with redundant components becomes: Aircraft will be killed if pilot is killed or fuel tank is killed
or both the engines are killed i.e., kill equation is,
(PILOT) OR (FUEL TANK) OR [ENGINE1) AND (ENGINE 2)]
The equation (3.9) for the probability of aircraft survival given a random hit on it in this case is modified
to: Aircraft will be survive if pilot survives and fuel tank survive and one of the engines survives
i.e., probability of survival of aircraft is,
(PILOT).AND.(FUELTANK).AND.[(ENGINE1).OR.(ENGINE2)]
That is
PS/H = (P S/H p
.P
S /H f (1– Pk / H e1 . Pk / H e 2 ) ) ...(3.13)
78 System Modeling and Simulation
If we assume that single hit cannot kill both the engines, then all the components killed are mutually
exclusive. Thus from (3.13) one gets:
is modified as
PS/ho = Ps1(1 – Ps2 . Ps3)...Ps4...Psc ...(3.14a)
As earlier, converting survivability into vulnerability, one can get an equation similar to (3.11).
In this chapter only the basic technique of mathematically evaluating the various cases of aircraft
orientation has been discussed. All these cases can occur in a single aircraft when it manoeuvres.
We will use these results in sixth chapter and generalise the vulnerability mode.
12345678
12345678
12345678
12345678
12345678
4
12345678
12345678
12345678
DISCRETE SIMULATION
weapon system, techniques of pure mathematics are not John von Neumann (Neumann János) (December
so handy and one has to opt for other techniques. In 28, 1903 – February 8, 1957) was a Hungarian
mathematician and polymath of Jewish ancestry
this chapter we will be discussing a different technique who made important contributions in quantum
which is also quite versatile in solving various problems physics, functional analysis, set theory, economics,
computer science, numerical analysis, hydrody-
where events are random. This technique is called namics (of explosions), statistics and many other
Monte Carlo simulation. Computer simulation is one of mathematical fields.
the most powerful techniques to study various types of Most notably, von Neumann was a pioneer of the
modern digital computer and the application of
problems in system analysis. The conceptual model is operator theory to quantum mechanics (see Von
the result of the data gathering efforts and is a Neumann algebra), a member of the Manhattan
Project Team, and creator of game theory and the
formulation in one’s mind (supplemented by notes and concept of cellular automata. Along with Edward
diagrams) of how a particular system operates. Building Teller and Stamslaw Ulam, von Neumann worked out
key steps in the nuclear physics involved in
a simulation model means this conceptual model is thermonuclear reactions and the hydrogen bomb.
converted to a computer model (simulation model).
Making this translation requires two important transi-
tions in one’s thinking. First, the modeller must be able to think of the system in terms of modelling
paradigm supported by the particular modeling software that is being used, second, the different
possible ways to model the system should be evaluated to determine the most efficient yet effective
way to represent the system. Many problems which can not be solved by mathematical methods can
be solved by Monte Carlo simulation.
80 System Modeling and Simulation
The term Monte Carlo was introduced by John von Neumann and Stamslaw Ulam during World
War II, as a code word for a secret work at Los Alamos. Monte Carlo is the name of a city in Monaco
which is famous for gambling casinos. You know that gambling is basically based on random numbers
or throwing of dice, which is also another form of random number generator. Problems dealing with
dice were discussed in second chapter. In Monte Carlo technique random numbers are generated and
used to simulate random events in the models. Aim of the present chapter will be to discuss various
methods of random number generation and apply them to physical problems. Later on, some case studies
such as ground target damage, will also be taken up.
Missile has proved its worth in destroying enemy targets with sufficient effectiveness involving
lesser cost, as compared to that by aircraft attack. Aircraft, although highly accurate are quite costly
and involve enemy attrition rate (i.e., damage of aircraft due to enemys fire power), whereas missile
has a capability of delivering sophisticated heavy warhead loads, with sufficient accuracy, without
much attrition rate. Apart from this, risk of human life loss is also involved in the attack by aircraft,
which is not present in the case of missile attack (see chapter ten on Cost effectiveness studies).
Some of the simple models of damage due to missile attack are discussed in this chapter using Monte
Carlo simulation. Aim of the present study is to demonstrate the power of simulation technique and
its applications to various practical problems. Case studies like runway and battle field denial by missile
warheads will be taken up.
In section 4.1, general introduction and the method of generation of random numbers is discussed
alongwith few case studies. Apart from studying the methods of generation of random numbers it
is important to study whether these numbers are really random. In section 4.2, methods to test the
randomness of numbers are discussed. In section 4.3, an inverse transform method for the generation
of random numbers will be studied and in section 4.4, different methods of generating normal random
numbers are discussed and few simple case studies will be taken up. In sections 4.5, and 4.6, of
this chapter, Monte Carlo random number generation method will be used to simulate the missile attack
and to assess the damage incurred to a runway and a battle field.
Generally question is raised, how authentic is the simulation technique? Because of various
methods available for the generation of random numbers, one always has a doubt about the
authenticity of the method. Each method will have its own errors. Then how to validate the model?
Even mathematical model will involve lot of assumptions, which may distract us from ground
realities. Another way is to compare the results with available experimental data. All these issues
will be studied in this chapter.
square. But these points are not random. By fixing the coordinates of one of the corner of square,
we can determine the coordinates of all these points. These points are uniformly distributed but
not randomly. Now we generate ten points in the bigger square, by some random number generation
technique. If we observe that approximately one point out of these lie in each small square, then
in this case we will say points are uniformly distributed and are random, because their location
is not fixed and is unknown. When we get into the detailed study of uniform random numbers,
we will come to know, it is almost impossible to generate true uniform random points in an area.
In the following paragraphs, we will first try to learn, how to generate uniform random numbers
between two given numbers. This is one dimensional case and called Monte Carlo method of
generating random numbers.
Because the sampling from a particular distribution involves the use of uniform random numbers,
stochastic simulation is often called Monte Carlo simulation. Uniform random numbers are the
independent random numbers, uniformly distributed over an interval [0,1], and are generally available
as a built-in function in most of the computers.
Length = 10 units
Height = 4 units
After using Monte Carlo simulation to test 10,000 random points, we will have a pretty good
average of how often the randomly selected location falls within the black area. We also know from
basic mathematics that the area of the rectangle is 40 square units [length × height]. Thus, the black
area can now be calculated by:
number of black hits
Black area = × 40 square units
10, 000 hits
Naturally, the number of random hits used by the Monte Carlo simulation doesn’t have to be 10,000.
If more points are used, an even more accurate approximation for the black area can be obtained.
#include<stdlib.h>
#include <math.h>
main( )
{
int a,b,k,j,i,m,nn,seed,r[50];
cout<<”\nEnter the integer value of a,b,m, leave space between a,b,m”;
cin>>a>>b>>m;
cout<<”\nEnter the integer value of seed”;
cin>>seed;
cout<<”\nEnter the number of Random numbers to be generated”;
cin>>nn;
r[0]=seed;
for(i=0; i<=nn; ++i)
{
r[i]=(a*r[i-1]+b)%m;
cout<<r[i];
}
/* Program ends*/
}
{
long int i,s,z1,x,nd,seed;
int n;
float z,y;
cout<<“\n Give the seed number of four digits”;
cin>>seed;
cout<<“\n Give the number of random numbers to be generated”;
cin>>n;
for(i=0; i<=n; ++i)
{
y=int(seed*seed/100.0);
//cout<<y<<endl;
z=(y/10000.0);
z1=int(z);
x=int((z–z1)*10000);
seed=x;
cout<<x<<“\n”;
}
}
OUTPUT OF PROGRAM
A sequence of 25 random numbers with seed value 3459 is given below.
9645 0260 0675 4555 7479 9353 4785 8962 3173
678 4596 1231 5153 5534 0250 0625 3905 2489
1950 8025 4005 0400 1599 5567 9913 2675
#include <time.h>
main( )
{
int x=0, y=0, t=0;
//Parameters have been initialised at definition level.
while(y<10)
{rand( );
if(rand( )>=0 ||rand( )<=4)
y++;
else if(rand( )=5 )
y--;
else if(rand( )>=6 ||rand( )<=7)
x++;
else
x--;
}t++;
cout<<“t=”<<t<<“x=”<<x<<“y=”<<y;
}}
1 6 R 1 0
2 1 F 1 1
3 0 F 1 2
4 1 F 1 3
5 7 R 2 3
6 0 F 2 4
7 6 R 3 4
8 6 R 4 4
9 7 R 5 4
10 4 F 5 3
11 7 R 6 3
12 2 F 6 4
13 1 F 6 5
14 5 B 5 5
15 3 F 5 6
Discrete Simulation 87
It is assumed that he takes one step per minute and his destination is ten steps in the direction
of y-axis. Find out the time taken by the drunkard to reach the destination. In order to solve this
problem, we generate ten uniform random numbers, between 0 and 9 using Congruential method.
If number lies between 0 and 4, y-coordinate is incremented by one step. If number is 5, y-coordinate
is decremented by one step. If number is 6 or 7, x-coordinate is incremented by one step and if number
is 8 or 9, x-coordinate is decremented by one step. Algorithm of this logic is given in C++ program
4.3 and results are given in Fig. 4.1.
f(x) fmax
(p,q)
O
a
x b
Fig. 4.2: Rejection method.
88 System Modeling and Simulation
3. If q > f (y), then reject the pair (u, v), otherwise accept u as the random number following
the distribution f (x).
This method works as all the random numbers accepted lie below the curve y = f (x). Only
restriction is that this method works for random number in a limited area, which is bounded by lower
and upper limits. In this method one has to generate large number of uniform random numbers which
may take more time in generating the desired random numbers.
In case curve f (x) does not have a maxima viz., curve is concave upward, then above step no.2
becomes q = f (b) . v, rest of the procedure is same (Fig. 4.3).
In case there are more than one cusp in the probability distribution function, highest of maximum
values can be taken as fmax (see Fig. 4.4).
f(t)
f max
a t b
Fig. 4.4: Probability density function with two maxima.
Discrete Simulation 89
Let us denote by D, the maximum of D+ and D –. The critical value of D from Kolmogrov-
Smirnov tables for α = 0.5, and N = 10 is 0.410. Value of D from the Table 4.2 is max(0.14, 0.26)
i.e., 0.26 which is less than 0.41. Hence these random numbers are uniform with level of significance
α = 0.5.
90 System Modeling and Simulation
Another popular test for testing the uniformity level of random numbers is known as chi-square (χ2)
test. The chi-square distribution is a special form of Gamma distribution. The chi-square test provides
a useful test of goodness of fit, that is, how well data from an empirical distribution of n observations
conform to the model of random sampling from a particular theoretical distribution. If there were
only two categories say, success and failure, the model of independent trials with probability p of
success is tested using the normal approximation to the binomial distribution. But for data in several
categories, the problem is how to combine the test for different categories in a reasonable way. This
problem was solved as follows, by the statistician Karl Pearson (1857–1936). For a finite number
of categories m, let Ni denotes the number of results in category i. Under the hypothesis that the
Ni are counting results of independent trials with probability pi, for large enough n the so-called
chi-square statistics
( Ni – npi )
2
m
∑
i =1 npi
that is the sum over categories of (observed – expected)2/expected, has distribution that is approxi-
mately chi-square with m – 1 degree of freedom. In statistical jargon, a value of the statistic higher
than the 95th percentile point on the chi-square distribution with m – 1 degree of freedom would “reject
the hypothesis at the 5% level”.
Thus the chi-square test uses the sample statistic
(Oi − Ei ) 2
n
χ = ∑
2
...(4.4)
i =1 Ei
where Oi = Ni is the observed random numbers in i-th class, Ei is the expected number in the i-th
class and n is the number of classes. For the uniform distribution Ei, the expected number in each class
is given by
N
Ei =
n
for equally spaced n classes, where N is the total number of observations.
2
It can be shown that the sampling distribution of χ is approximately the chi-square distribution
with (n – 1) degree of freedom. To make the process clearer, let us generate hundred rn numbers of
uniform random numbers lying between 0 and 1. Divide these numbers in ten classes (n) of equal interval
so that random numbers less than or equal to 0.1 fall in the 1st class, those of 0.1 < rn ≤ .2 fall in
2nd class, 0.2 < rn ≤ .3 fall in 3rd class and so on. In this way, let us assume Oi number of random
numbers fall in the i-th class. Ei here is 100/10 = 10. Table 4.2 gives Oi in different classes in second
column. Third column gives the difference Oi – Ei and fourth column square of Oi – Ei. Then using
equation (4.4) we get,
χ2 = 32/10 = 3.2
2
From the χ2 tables in the Appendix 4.3, we find that for degree of freedom 9, value of χ
for 95% level of confidence is 16.919, which is more than our value 3.2. Thus random numbers
of Table 4.2 are uniform with 95% level of confidence.
Number of outcomes in column two of the above table are taken as expected values Ei as in
chi-square test. Exact outcome is taken as Oi. Like chi-square test standard deviation is computed.
Example 4.3: Poker Test
A sequence of 10000 five digits numbers are generated. Outcome of different categories of combination
is as given in Table 4.5.
The degrees of freedom in this case are six, one less than the number of combinations.
The critical value of χ2 for six degree of freedom at a = 0.01 is 16.8. The value obtained in
the last column is 12.8336, which is less than the critical value. Hence, the number generated
are independent.
Discrete Simulation 93
.33
.67
1.0
Example 4.4: Following random numbers have been generated by congruential method. Test their
randomness.
49 95 82 19 41 31 12 53 62 40 87 83 26 01 91 55 38 75 90 35 71 57 27 85
52 08 35 57 88 38 77 86 29 18 09 96 58 22 08 93 85 45 79 68 20 11 78 93
21 13 06 32 63 79 54 67 35 18 81 40 62 13 76 74 76 45 29 36 80 78 95 25 52.
Solution: These are 73 numbers and we divide them into 72 pair such that 49, 95, is first pair,
95, 82 is second pair, 82, 19 is third pair and so on. If we call each pair as R1 and R2 then condition
is R1 ≤ .67 & R2 ≤ 1.0 thus it goes to cell C(2, 3). Similarly we place other pairs in their respective
cells. Following table gives the frequency of pairs in each cell.
94 System Modeling and Simulation
Therefore
χ2 = 24/8 = 3.0
In this case since there are two variable R1 and R2 and hence the degree of freedom is nine minus
two i.e., seven. The criterion value of χ2 for seven degree of freedom at 95% level of confidence
is 14.067. This value is much higher than the values obtained by present test and hence given random
numbers are not auto correlated.
F ( x) = ∫ f (x)dx
−∞
...(4.5)
where 0 ≤ F ( x) ≤ 1, since the integral of probability function f (x) over values of x varying from
−∞ ≤ f ( x) ≤ + ∞ is 1. You can observe some similarity between random numbers generated in previous
section and CDF. Can you guess what? Both vary between 0 and 1. Figure 4.5 shows a typical cumulative
distribution function.
The cumulative function can be solved for x; i.e., if y = F (x), then
x = F −1 ( y ) ...(4.6)
Now we will prove that if y is uniformly distributed over a region 0 ≤ y ≤ 1 then variate X, has
a distribution whose values x, are given by equation (4.5).
Discrete Simulation 95
If X is a random variate governed by the given probability density function f (x), and that
Y = F(x) is a corresponding value of the cumulative distribution, then
P{Y ≤ y0} = P{X ≤ x0}
But by the definition of cumulative density function, P{ X ≤ x0 } = F{x0 } = y0, hence
P{Y ≤ y0 } = y0 ...(4.7)
which is the expression for the cumulative uniform distribution within the interval (0,1). Thus we
see that Y is uniformly distributed in the interval (0,1), irrespective of the distribution of X. This fact
is also clearly shown in Fig. 4.6.
Therefore in order to generate a random variate X with distribution f (x), we first generate a uniform
random variate U lying between (0,1). Then we obtain X as,
X = F −1 (U )
It is often observed that it is not possible to integrate the function f (x) to get the CDF (normal
distribution function is one of the example). Some time it is possible to get CDF of f (x) but is not
possible to invert it. Thus inverse transform method fails. In such cases we have to adopt alternative
methods. There are other methods for the generation of random numbers, for example, composition
method and acceptance-rejection method. Since study of basic theory of random number generation
is not the aim of the book, readers may see reference [53, 77, 86] for detailed study. A typical example
of determination of required random number is given below.
Example 4.5: Generate a random variable with uniform distribution f (x) given by,
1
, a≤ x≤b
F(x) = b – a
0, otherwise
Solution: Cumulative distribution function is given by,
0, x < a,
x–a ,
F(x) = a≤ x≤b
b – a
1, x>b
and thus solving for x = X, one gets,
X = F–1(U) = a + (b – a)U ...(4.8)
96 System Modeling and Simulation
Thus X is a uniform random variable lying between limits a and b. In the following example,
we illustrate a problem which is reverse of the problem given in example 4.5 i.e., when distribution
of random variate is given, determine the distribution function.
Example 4.6: Generate a distribution function f (x) when a random variable X as a function of
uniform distribution is given.
Solution : Let us consider a random variable given by X = –2ln (U), where U is a uniform random
variable between [0,1]. Our aim is to find distribution function of variable X.
Now by definition (equation 4.3), cumulative distribution function of random variable X is given
by
FX (x) = P( X ≤ x)
or
P( X ≤ x) = P( X ≤ −2ln(U ) ≤ x) = P(2 ln(U ) ≥ − x) = P(ln(U ) ≥ − x / 2)
= P (eln(U ) ≥ e − x / 2 )
∞
∫ fU (u ) . du
–x/2
= P (U ≤ e )=
e– x / 2
where fU (u) is the probability density function of a random variable that is uniform on [0,1]
x
1 –( z – µ )2 / 2 σ2
FX (x) =
2πσ ∫e
–∞
dz
We know that there is no closed-form expression for FX (x), and also it is not possible to solve
this expression for random variable X.
One of the approximation for FX (x) is given by Shah (see reference [18]), which states that,
for µ = 0 and σ = 1,
FX (x) = x(4.4 – x )/10 + 0.5 ...(4.9)
with an error of the order 0.005 for x lying between 0 and 2.2. This method is however not
recommended for the generation of normal random numbers.
Algorithm that generates N (0,1) random variables can always be modified to become general
N(µ, σ 2) generators, as follows:
Y=
∑ i=1Ui − k/2 ...(4.12)
k /12
has mean 0 and variance 1. The distribution function of Y is close to that of the standard normal
distribution for large k. In fact in equation (4.11), k need not be very large. According to Stirling’s
approximation k = 12 is sufficient to generate normal number accurately.
98 System Modeling and Simulation
X1 = −2ln(U1 ) sin(2πU 2 ),
is exactly N(0,1).
Using this method, one can generate a sequence of independent N( µ, σ 2) random variables to be
This method is relatively fast, and provides exact normal random variables. Sometimes cos (.)
in place of sin (.) is also used. This also gives normal random numbers different from those given
by sin (.) and is independent of each other. Care must be taken that all Us’ must be independently
chosen from different streams of independent numbers.
1
fmax =
2πσ 2
equal to the lethal radius of the bomb. Then, another stream of pair of random numbers, from an
appropriate uniform distribution, are generated representing target points on the target. If these points
lie in one of the three circles, they are counted, else rejected. Let M of these lie with in one or more
of the lethal circles with their centers at simulated points of impact. Also, let K represents the number
of points, out of M, which lie within at least two lethal circles. Then, the expected coverage of the
target can be estimated by the sample mean E of
E = M/N
and the expected overlapping coverage, by the sample mean D of
D = K/M
Using relevant techniques of statistical inference, as described below, we determine the error in
the estimated coverage by Monte Carlo Simulation method.
Let j simulation runs be made using different streams of random numbers for simulating the given
warheads impact points and N target-points. The sample means E may be considered to estimate
the expected target coverage and the associated accuracy may be specified by the confidence interval
with 100(1– α)% level of confidence. Using technique of statistical inference [79], the end points
of this confidence interval turn out to be
α σ
E ± t , j –1 . E ...(4.16)
2 j
where t(α/2, j – 1) represents the 100(1 – α)% percentage point of the Student’s t-distribution
with (j – 1) degrees of freedom; and σ E2 , the sample variance. Thus, when the (true) expected
coverage is estimated by E, the maximum error δ of the estimate with confidence 100 (1– α)%
would be given by (see appendix 4.2)
σE
δ = t ( α /2, j − 1) . ...(4.17)
j
It may be noted that j represents the number of observations (simulation runs) and must be
distinguished from N which may be regarded as the length of the simulation run. Similarly, the
overlapping coverage D may be analyses.
will extend this model, and study the criteria of airfield denial. Monte Carlo technique of simulation
is used to find the number of missiles required to be dropped on the runway tracks given in the
section 1.3.1, to ascertain a specified level of damage.
Damage to airfield is to such an extent that no aircraft can land or take off from it. It is well
known that modern aircraft can land or take off even if a strip of dimensions 15 × 1000m is available.
In order to achieve this we choose few aim points on the runway. These aim points are called Desired
Mean Point of Impacts (DMPI). These aim points are taken as the centre of the Desired Mean Area
of Impacts (DMAIs), which are demarcated on the runway. Let (xd , yd ) be the co-ordinates of one
of the aim points. To find impact point, two normal random numbers x and y are generated (using
Box-Mullar method)
x = −2 log(u1 ) sin (2πu2)
bx – x g + by – y g
i 1
2
i 1
2
≤ ( Rwh – rb )
RH DENIED -1
TTDENIED - 1
ABH DENIED - 1
NTALAS - 10
NSUODESS - 6
Fig. 4.9: A computer output of runway denial model using new denial criteria.
102 System Modeling and Simulation
Knowing the position of all the bomblets, it is ascertained that each strip of width Ws has at least
one bomblet. If all the DMAIs, are denied, experiment is success otherwise failure. Trial is repeated
for large number of times (say 1000 times) and the probability of denial is calculated as the ratio
of the number of successes to the number of trials. To ascertain the correct probability of denial
programme has been run n times (say 15 times) and the actual probability of denial has been obtained
as the average of these n probabilities. Figure 4.9 gives the output of the compute model.
EXERCISE
APPENDIX 4.1
Y = –2 ln (U1 ) cos(2 p U 2 )
104 System Modeling and Simulation
APPENDIX 4.2
σx =
∑ ( x – µ x )2
N
where N = total number of possible samples. But it is not possible to take total number of possible
samples for determining s x , if population is too large. Thus limited number of samples of finite size
are chosen out of the large population, and standard error is estimated as
σ N–n
σx =
n N –1
where σ = the population standard deviation,
N = population size,
n = sample size.
If population is infinite, standard deviation can be calculated as
σ
σx =
n
When the population standard deviation (σ) is known, we can directly compute its standard error
of the mean. Thus, the interval estimate may be constructed as
x – zσ x < µ < x + zσ x
x−µ
where z = under the normal curve. For example, for 95% confidence, coefficient z = 1.96.
σ
Discrete Simulation 105
In many situations, not only is the population mean unknown but also the population standard
deviation is unknown. For such a case, it appears intuitively that the sample standard deviation (s)
is an estimation of the population standard deviation to σ. The computation of s and σ are given by
∑ (x − x ) 2
∑ ( x − µ) 2
s = , σ =
n N
But although there is a similarity between s and σ, we must remember that one of the important
criteria for a statistic to qualify as an estimation is the criterion of unbiasedness. The sample standard
deviation is not an unbiased estimator of the population standard deviation.
The unbiased estimator of the population standard deviation is denoted by s and is given by
∑ (x – x ) 2
σ̂ = n –1
Thus, if the population standard deviation is unknown, the sampling distribution of means can
be assumed to be a approximately normal only when the sample size is relatively large (> 30). In
this case interval estimate for large samples is altered slightly and is written as
x – z σˆ ∠ µ ∠ x + z σˆ x
x – t (α / 2, n – 1) σ x ∠ µ ∠ x + t (α / 2, n – 1) σ x
Here α is the chance of error and n is degree of freedom in t-distribution (sample size). Like
the z-value, the value of t depends on the confidence level. For t-distribution see reference to some
standard text book on statistics [54,79].
106 System Modeling and Simulation
APPENDIX 4.3
APPENDIX 4.4
Degree of
D0.10 D0.05 D0.01
freedom (N)
APPENDIX 4.5
INDICATOR FUNCTION
In mathematics, an indicator function or a characteristic function is a function defined on a set
X that indicates membership of an element in a subset A of X.
The indicator function of a subset A of a set X is a function
IA : X →{0,1}
defined as
1 if x ∈ A,
IA(x) =
0 if x ∈ / A.
12345678
12345678
12345678
12345678
12345678
5
12345678
12345678
CONTINUOUS SYSTEM
SIMULATION
So far we have discussed problems and techniques which are stochastic in nature. In the present chapter,
problems of continuous nature will be studied. Modeling and simulation of various problems like pursuit
evasion of aircraft, fluid flow, flight dynamics and so on come under continuous system simulation.
Continuous system generally varies with time and is dynamic. As it has been mentioned earlier too, there
is no fixed method for modeling a particular problem. However, one of the basic tools for modeling
beyond doubt is mathematics. One can model any situation with the help of mathematics. Sometime
mathematical models cannot be solved analytically, as we have seen in the case of hanging wheel of
a vehicle. There was a time when approximations were made to simplify the mathematical equation of
the problem. But now with the advances made in modern computers, any continuous model can be
converted to digitization and worked out with the help of computer programming. We had assumed
in this book that reader is conversant with computer programming with special reference to C++ language.
In first few sections, we will discuss the basic mathematical tools to be used in continuous
simulation. It is not possible to study these tools in details, as full book is needed for each topic,
but attempts will be made to cover the topic briefly and reference to required book will be made
wherever required.
∂
dm = (ρdx . dy . dz) .dt
∂t
∂ρ . . .
= dx dy dz dt ...(5.1)
∂t
since x, y, z, and t are independent coordinates in three-dimensional space.
dz
dy
dx
Change in the mass of fluid in the cube is due to the movement of fluid in and out of its surfaces.
Let velocities of the fluid in x, y, z directions be u, v and w respectively, where these velocity
components are function of space coordinators x, y, z and time t. If in time dt, mass (ρu )x enters
from the surface dy . dz and (ρu )x + dx leaves from the other parallel face then,
∂
[(ρu ) x − (ρu ) x + dx ] dt dy dz = – (ρu) dt dxdy dz ...(5.2)
∂x
where higher order terms in the expansion of (ρu ) x + dx have been neglected, being comparatively small.
Similarly we get two more relations for movement of fluid along y- and z- directions i.e.,
∂
[(ρv) y − (ρv) y + dy ] dt dxdz = − (ρv) dt dx dy dz
∂y
∂
(ρw) dt dx dy dz
[(ρw) z − (ρw) z + dz ] dt dx dy = – ...(5.3)
∂z
Sum of the fluid motion out of three surfaces of the cube is equal to rate of change of mass
in side the cube. Thus from equations (5.1), (5.2) and (5.3) one gets,
∂ρ ∂ ∂ ∂
+ (ρu ) + (ρv) + (ρw) = 0 ...(5.4)
∂t ∂x ∂y ∂z
This equation is the equation of continuity. This result in vector form can be written as,
∂ρ
+ div(ρv ) = 0
∂t
where the velocity v has the components u,v, w in cartesian coordinates and div is the divergence.
112 System Modeling and Simulation
∂u ∂u ∂u ∂u ∂P
ρ + ρu + ρv + ρw =− ...(5.5a)
∂t ∂x ∂y ∂z ∂x
Similarly two equations in y- and z-directions are obtained as,
∂v ∂v ∂v ∂v ∂P
ρ + ρu + ρv + ρw =− ...(5.5b)
∂t ∂x ∂y ∂z ∂y
∂w ∂w ∂w ∂w ∂P
ρ + ρu + ρv + ρw = – ...(5.5c)
∂t ∂x ∂y ∂z ∂z
These equations when written in vector form are,
dv ∂v
ρ = ρ + ρ(v .grad)v = – grad P ...(5.6)
dt ∂t
d 1
ρ E + (u 2 + v 2 + w2 ) dt dx dy dz
dt 2
Where the total time derivative is used to account for the displacement of the element during
the interval. This change in energy must be equal to the work done by the fluid on the surfaces
Continuous System Simulation 113
of element volume. The work done in time dt on an area dydz in motion along x-direction, is the
product of force and displacement, or P u dtdydz, and net amount of work done on the two faces
of the volume is
∂( Pu )
[( Pu ) x − ( Pu ) x+ dx ] dtdydz = −
dt dx dy dz
∂x
The work done on other faces is obtained in the same way, and equating the total work on all
the faces to the change in energy, one gets
d 1
E + (u 2 + v 2 + w 2 ) = − ∂ ( Pu ) + ∂ ( Pv) + ∂ ( Pw)
ρ ...(5.7)
dt 2 ∂x ∂y ∂z
which is the equation for energy conservation. This equation in vector notation is given as,
d 1
ρ E + (v . v ) = – div( Pv) ...(5.7a)
dt 2
where v = (u , v, w).
The equations (5.4),(5.6) and (5.7) are non-linear equations and can not be solved by analytic
methods. Therefore, one has to go for numerical techniques. Some of the techniques for numerical
computations and their programs in C++ language are given in appendix 5.1 for the convenience of
reader. Those who want to study this subject in details can refer to some book on computational
fluid dynamics.
d2x D dx K K
2
+ + x = F (t )
dt M dt M M
which can be written as
d2x dx 2
+ 2ζω + ω 2 x = ω F (t ) ...(5.8a)
dt 2 dt
where 2ζω = D / M, ω2 = K / M
114 System Modeling and Simulation
We integrate this equation using Runge-Kutta method of fourth order (see appendix 5.1). In
order to use this method, equation (5.8a) is to be converted in first order equations. Equation (5.8a)
is converted to two first order differential equation as follows,
dx
=y
dt
dy
= f ( x, y, t ) = ω2 F (t ) – 2ζωy – ω2 x ...(5.8b)
dt
These are two homogeneous first order differential equations in x and y. Runge-Kutta method
gives us
m1 = yi
m1h
m2 = yi +
2
mh
m3 = yi + 2
2
m4 = yi + m3h
n1 = f ( xi , yi )
n1h mh
n2 = f ( xi + , yi + 1 )
2 2
nh m2h
n3 = f ( xi + 2 , yi + )
2 2
n4 = f ( xi + n3 h, yi + m3 h)
m + 2 m2 + 2 m3 + m4
xi +1 = xi + 1 h
6
n + 2 n2 + 2 n3 + n4
yi +1 = yi + 1 h ...(5.9)
6
h = ∆t
Using the above algorithm, and assigning some values to F(t), D, and K one can compute the
results of equation (5.8b). Figure 5.2 gives the variation of x verses ωt, when a steady force F is
applied to wheel at time t = 0.
Figure 5.2 shows how x varies in response to a steady force applied at time t = 0. Solutions
are shown for different values of ζ. It can be seen from the curves for different values of ζ that
oscillations of the wheel increase as ζ decreases. For oscillations to be minimum, ζ ≥ 1, that means
D2 ≥ 4MK. Here ζ is called the damping force. When motion is oscillatory, the frequency of oscillation
is determined from the formula
ω = 2π f
where f is the number of cycles per second.
Continuous System Simulation 115
Table 5.1: Lethal radii computed from the Scaling Law for various types of damage
Targets and Persons Persons Soft skinned Window glass Antitank Antipersonnel
types of (kill) (ear drum vehicles (breaking) mines mine
damage burst) (damage)
Pressure
required 4.5332 2.0332 1.3132 1.1332 7.6632 1.6632
(kg/cm2)
One of the most important property of shock waves is, the sudden increase in pressure behind it’s
front. This pressure (called peak over pressure) is mainly responsible for damage to structures. Some
of the important parameters of shock wave are: peak over pressure (P1 – P0), time duration of positive
pressure τ, time duration of negative pressure τ' and shock impulse I. Here P1 is the pressure behind
the shock front and P0 is ahead of it. These parameters are responsible for the determination of damage
due to blast. In Table 5.1, we give pressure range which can cause various damages.
zb
t
I = g
P1 – P0 dt ...(5.10)
0
where ( P1 – P0 ), is a function of time t and is of the form P1 = Pm e – at. Impulse is an important parameter
for the damage.
three circles in 5.4a and 5.4b are the positions of sound waves starting respectively from points A,
B and C at time when aircraft reaches point O. For details see Singh (1988). If we draw a tangent
to all these circles, it will pass through point O. Thus envelope of these circles is a cone with half
angle α. In Fig. 5.5, it can be seen that aircraft moves with speed v which is greater than sound
speed a, thus all the disturbance due to the motion of the aircraft remains confined in a cone whose
semi-vertex angle α is given by,
α = sin–1(1/M)
Fig. 5.4: Disturbance due to a body moving with sonic and supersonic velocity.
where M= v/a. Ratio M is called Mach number. Derivation of this equation is clear from Fig. 5.4b.
By the time aircraft moves from point A (B or C) to O, sound wave produced at A (B or C) reaches
the point A′ (B′ or C′), where line OA′ is tangent to all the three circles with centers at A,B and
C respectively. Thus,
A A ′ at 1
sin α = = =
OA ′ vt M
Ouch
Don’t hear
any noise
Fig. 5.5: Aircraft disturbances cannot be heard outside the Mach cone.
(Courtesy G. Gamow; frontiers of Physics)
118 System Modeling and Simulation
Thus we see that disturbance in case of supersonic motion remain confined in side the Mach
cone with semi-vertex angle α (see Fig. 5.5). It is known that when a supersonic aircraft is seen
by a person on the ground, its sound is not heard, but when it has crossed the person’s head, suddenly
a roaring sound is heard. This is because earlier, person was out of the range of the Mach cone and
heard no noise but as soon as aircraft crossed over his head, he came inside the Mach cone, in which
whole disturbance is confined and heard the roaring sound. Surface of the Mach cone is nothing but
the shock front. As we have stated earlier, the pressure behind the shock front is very high, by which
even buildings are damaged. This shock wave is weaker than the blast wave (that produced by
explosives), that is the reason it has been given different name after a famous scientist Ernest Mach
and is called Mach wave.
where
yb ( t ) – y f ( t )
tan( φ) =
x b (t ) – x f (t )
This problem although simple can not be solved analytically. We adopt the technique of numerical
computation. Conditions are, when distance (dist) ≤ 10km, missile is fired and bomber is destroyed.
If this distance is not achieved in 12 minutes from the time of detection, bomber escapes. We have
assumed that velocity of bomber is 20 km/m and that of fighter is 25 km/m. It is observed that target
is killed in six minutes with given velocities. Results of computation are given in Fig. 5.6 and a computer
program is given as program 5.1.
theta = atan((yf–yb)/(xb–xf));
cout<<“t=” << t <<“\n”;
cout<<“xb=”<<xb<<“yb=”<<yb<<“xf=”<<xf<<“yf=”<<yf<<“\n”;
cout<< ‘\t’<<“dist=”<< dist <<‘\t’<<“theta=” <<theta<<endl;
cout<<“\n”;
if ( t> tlim) {cout << “target escapes, game over”;break;}
else if(dist<=10.00)
{cout<<“target killed\n”;break;}
else
continue;
}
}
I
θ0 + D θ 0 + K θ0 = K θi ...(5.13)
Dividing this equation by I on both sides and making the following substitutions,
D K
2ζω = , ω2 = , ...(5.14)
I I
Thus equation (5.13) reduces to,
+ 2ζωθ + ω 2θ = ω 2θ
θ ...(5.15)
0 0 0 i
This is the second order differential equation and is same as equation (5.8a), and its solution is
provided in section (5.3). The result shows that the aircraft will have oscillatory motion for ζ ≤ 1 ,
which implies,
D 2 ≤ 4KI
Continuous System Simulation 121
There can be hundreds of example in continuous simulation and it is not possible to discuss all
of them. We have studied some of these, just to illustrate, what is continuous simulation.
d2x 1
m 2 =
– ρSCDV 2 cos(θ) ...(5.16)
dt 2
d2 y 1
m 2 =
− ρSCDV 2 sin(θ) + mg ...(5.17)
dt 2
θ = tan –1 ( dy / dx ) ...(5.18)
where
θ = angle of elevation,
CD = drag coefficient,
ρ = density of air,
g = acceleration due to gravity,
m = mass of the body.
It is known that drag coefficient is a function of the velocity of the projectile.
122 System Modeling and Simulation
EXERCISE
APPENDIX 5.1
Y-axis
yi
(xi, yi ) (xi+1, yi )
xi xi+1 X-axis
In the Fig. A.1, point xi+1 is nothing but x + ∆x, and xi,yi+1 is xi, yi + ∆ y. Using this nomenclature,
we can express differential equation (A.1) as
dy ∆ y yi +1 − yi
= =
dx ∆ x xi +1 − xi
Thus differential equation (A.1) becomes at i-th grid as
yi +1 = yi + ∆ x(axi + byi + c)
If initial values are given as when x = x0, y = y0, then
y1 = y0 + ∆ x(ax0 + by0 + c)
y 2 = y1 + ∆ x ( ax1 + by1 + c )
Y
P
Q
T
P
xi xi + 1 X
Fig. A.2: Errors in Euler’s method.
Now x1 = x0 + ∆ x, one gets y2 by substituting the values of y1 and x1 in equation for y2. Same
way we can compute the values for y3, y4,… . This method is the simplest one and is known Euler’s
method. We can easily see that in this method each successive value of dependent variable depends
on the previous value. Hence inaccuracies cropped in yi are propagated to yi+1 and slowly and slowly,
it deviates from the actual solution. In Fig. A.2, point T is (xi, yi) and point P is (xi + 1, yi) where
as point Q is (xi + 1, yi+1). But value y at (i + 1)-th point of grid, calculated by Euler’s method (tangent
at point T ) is at P′ which is less than the actual value of yi +1. Thus the error in yi+1 that is QP′
is added to each step and ultimately curve computed by Euler’s method deviates from actual curve.
In section A.3 we give an improved method which gives better approximation as compared to Euler’s
method. A computer program of Euler’s method is given below.
{
int i,n;
float x,y,xp,h,dy;
float func(float, float);
cout<<“\nSolution by Euler’s Method\n\n”;
/* Reading Initial data*/
cout<<“\n Input Initial values of x and y\n”;
cin>>x>>y;
cout<<“input x for which y is required\n”;
cin>>xp;
cout<<“Input step size h\n”;
cin>>h;
/*Compute number of steps required*/
n=(int)((xp–x)/h+0.5);
/*Compute y recursively at each step*/
for(i=1;i<=n;i++)
{dy=h*func(x,y);
x=x+h;
y=y+dy;
cout<<i<,x<<y<<“\n”;
}
/*Write the final results*/
cout<<“Value of y at x=”<<x<<“is”<<y<<“\n”;
}//End of main()
float func(float x, float y)
{
float f;
f =2.0 *y/x;
return(f)
}
dny
All numerical techniques for solving differential equations of the type = f (x, y) involve a series
dxn
of estimation of y (x) starting from the given conditions. There are two basic approaches that could
be used to estimate the value of y (x). They are known as one-step methods and multiple step methods.
In one step methods, we use the information from only one preceding point i.e., to estimate the
value yi, we need the condition at the previous point yi– 1 only. Multi step methods use information
at two or more previous steps to estimate a value. In this section we will discuss some of the integration
techniques for ordinary differential equations.
where y n (x0) is the n-th derivative of y (x), evaluated at x = x0. The value of y (x) is known if we
know its derivatives. To understand the method let us consider a differential equation,
dy
= y′ = 2 x2 + 3 y3 + 5 (A.3)
dx
with the initial conditions y(0) = 0.5 at x = 0. Now
y′′ = 4 x + 9 y 2 y ′
y′′′ = 4 + 18 y( y' )2 + 9 y 2 y ′′
At x = 0, y (0) = 0.5 and therefore
y′(0) = 5.625,
y′′(0) = 9 (1) (8) = 12.656,
y′′′(0) = 4 + 18(0.5) (5.625)2 + 9.(0.25) (12.656)
= 4 + 284.765 + 28.476 = 317.241
Substituting these values in the Taylor’s expansion (A.2),
y(x) = 0.5 + 5.625x + 6.328x2 + 52.87x3 +…
Number of terms used in the above equation depend on the accuracy of y required. There are
some shortcomings in this method. It is not very accurate as sometimes large number of terms have
to be considered and for higher order terms, derivatives become more and more cumbersome.
yi +1 = yi + hf xi + xi +1 yi + yi +1
, , h = ∆x
2 2
= yi + hf ( xi + h / 2, yi + ∆y / 2) (A.4)
∆y is the estimated incremental value of y from yi and can be obtained using Euler’s formula as
∆y = hf ( xi , yi )
Thus equation (A.4) becomes
yi +1 = yi + hf ( xi + h / 2, yi + hf ( xi , yi ) / 2)
= yi + hf ( xi + h / 2, yi + hm1 / 2), (A.5)
Continuous System Simulation 127
where
m1 = f ( xi , yi )
= yi + m2 h
where
m2 = f ( xi + h /2, yi + hm1 /2)
this method is called Modified Euler’s method or improved polygon method. Sometimes it is also
called mid point method.
Let us integrate equation
dy
= 2 y/x
dx
with initial conditions at x = 1, y = 2.0
by Polygon method and compared with its integration with Euler’s method.
Now we assume that h = 0.25 then
y (1) = 2.0
{
int i,n;
float x,y,xp,h,m1,m2;
float func(float, float);
cout<<“\n Solution by polygon Method\n\n”;
/* Reading Initial data*/
cout<<“\n Input Initial values of x and y\n”;
cin>>x>>y;
cout<<“input x for which y is required\n”;
cin>>xp;
cout<<“Input step size h\n”;
cin>>h;
/*Compute number of steps required*/
n=(int)((xp–x)/h+0.5);
/*Compute y recursively at each step*/
for(i=1;i<=n;i++)
{
m1=func(x,y);
m2=func(x+0.5*h,y+0.5*h*m1);
x=x+h;
y=y+m2*h;
cout<<i<,x<<y<<“\n”;
}
/*Write the final results*/
cout<<“Value of y at x=”<<x<<“is”<<y<<“\n”;
}//End of main()
float func(float x, float y)
{
float f;
f =2.0 *y/x;
return(f)
}
h mh
m2 = f ( xi + , yi + 1 )
2 2
h mh
m3 = f ( xi + , yi + 2 )
2 2
m4 = f ( xi + h, yi + m3 h) (A.6)
m + 2m2 + 2m3 + m4
yi +1 = yi + h 1
6
These equations are called Runge-Kutta method of fourth order. For detailed derivation of the
method readers may refer to [4].
∂u ∂u
+ a (u , x, t ) = 0, 0 ≤ x ≤ 1, t > 0, a > 0 (A.7)
∂t ∂x
where initial values are
u ( x, 0) = F(x), 0 ≤ x ≤1
u (0, t ) = g(t), 0≤t ≤T
where f (x) represents the initial conditions and g(t) the boundary conditions at x = 0. No boundary
conditions are required at x = 1.
A difference method might be set-up by selecting a grid with spacing ∆ x and ∆ t in the x and t
direction respectively. Since in the present case there are two independent variables x and t, we use
u(x, t) as u nj = u( jDx, nDt ), where superscript is for n-th point along t-axis and subscript is for j-th point
along x-axis in the grid. There are many ways of converting the equation (A.7) into difference equation.
A popular choice for this simple transport equation is the following difference equation.
u nj +1 − u nj u nj − u nj −1
+ a(u nj , x j , t n ) = 0 (A.8)
∆t ∆x
which can be solved at all positive integer values of j and n once values have been given for
j = 0 and n = 0. As a special case if a is a constant, f (x) = sin(x) and g (t) = – sin(at), then
the solution to the equation (A.7) is u(x, t) = sin(x – at). The difference equation (A.7) can be
solved by letting u 0j = f ( j ∆x) and letting u0n = g (n ∆t ).
For detailed study of numerical techniques using difference equations, readers are referred to
reference [4].
130 System Modeling and Simulation
dny
All numerical techniques for solving differential equations of the type = f ( x, y) involve a series
dx n
of estimation of y (x) starting from the given conditions. There are two basic approaches that could
be used to estimate the value of y (x). They are known as one-step methods and multiple step methods.
In one step methods, we use the information from only one preceding point i.e., to estimate the
value yi, we need the condition at the previous point yi– 1 only. Multi step methods use information
at two or more previous steps to estimate a value. In this section, we will discuss some of the integration
techniques for ordinary differential equations.
12345678
12345678
12345678
12345678
12345678
6
12345678
12345678
12345678
One of the applications of dynamic simulation is in the field of aircraft vulnerability studies. This is
an important field of aircraft design industry. The study of vulnerability of a combat aircraft against
ground based air defence system is of utmost importance for the design and development of aircraft,
so that it is capable of surviving against it. Also such a software is used for the evaluation of various
anti aircraft weapon system which are under development. This study is also needed in planning the
Air Defence System (ADS) against an enemy attack. So far damage to stationary targets and
mathematics behind it had been discussed in the previous chapters. In chapter three, it was assumed
that projection of aircraft and its vital parts on a given plane have been provided as inputs to the model.
Also effect of motion of aircraft is also ignored. In this chapter study of single shot hit probability
to a moving target will be undertaken. For this purpose, aircraft will be taken as a moving target.
Preliminary study of vulnerability, effect of redundancy and overlapping has already been given in
chapter three. These concepts will be used in the present chapter. Number of models dealing with
vulnerability of aerial targets, have been reported in the literature [2, 38, 51, 66–69]. A simplistic form
of model where areas of vulnerable parts projected on a given plane are known, has been given by
Ball [5] (see chapter three). There it was assumed in that the projection of aircraft as well as its vital
parts on a typical surface is available as an input data. In the present chapter, we will give complete
algorithm for the evaluation of vulnerability of a typical aircraft. Only inputs used in this model are
the structural data of the aircraft and characteristics of weapon system.
A dynamic model of aircraft vulnerability due to ground air defence system, where aircraft is
assumed to be moving along an arbitrary profile, is considered. A proximity fuse shell, fired from
ground towards aircraft, has been assumed to explode in the near vicinity of the target aircraft. This
shell has fragments as well as blast effects. If it explodes near the target, main damage may be due
to blast only. Damage to aircraft body due to explosive charge as well as fragments, when warhead/
ammunition explodes in its near vicinity has been considered in this chapter. Kill criterion has been
taken as the minimum number of fragments required to penetrate and kill a particular part based on
132 System Modeling and Simulation
the total energy requirement. In the case of damage due to blast waves, it is assumed that the probability
of kill is one, since energy released by shell is sufficient for catastrophic explosion. This assumption
is based upon the total impulse transmitted to the aircraft structure. A typical aircraft and a typical
air defence gun with DA/VT fused ammunition has been considered for the validation of the model.
Structural data of aircraft as well as its vital parts is taken in the form of triangular elements [66].
These triangular elements are obtained by dividing surface of whole aircraft and it’s vital parts into
small three dimensional triangles. For this purpose one has to go to the drawings of the aircraft and
obtain the (x, y, z) co-ordinates of apexes of all the triangles. Kill criterion due to fragment hits has
been modified and is based on fragment energy concept. For DA fused impact, it has been assumed
that shell first hits the outer surface of aircraft and then explodes. Kill in this case is mainly due to
impact and explosive energy of the shell. A three dimensional model for single shot hit probabilities
has been presented in this chapter for the case of proximity fused ammunition. The effect of redundant
vulnerable parts on the overall kill probability of aircraft is also studied in this model.
criterion due to fragment hits is based on fragment energy concept. Fragment energy concept is
explained in section 6.5.1. Damage to aircraft body due to explosive charge as well as fragments,
when warhead/ammunition explodes in its near vicinity has been considered. Aircraft in fact is a
complex system, in which each and every part has an important role to play. It has large number
of vital components. But for the sake of simplicity only five vital parts are considered. These
components are avionics, pilot 1, pilot 2, two fuel tanks and two engines. Fuel tanks and engines
are taken as redundant vital parts. Fuel tanks are of two types; one is fuselage fuel tank and other
is in two wings. Two fuel tanks located in two wings are treated as one part. Two pilot are assumed
to be non-redundant parts.
Probability of kill of aircraft/vital part depends on probabilities of its detection, hit and fuse
functioning. Probability of kill of j-th vital part of an aircraft by a round is defined as [66].
j j
Pk j = Pd P h Pf P k / hf ...(6.1)
where Pd , Phj , Pf , P j k / hf respectively are probabilities of detection of the aircraft, hit on the j-th vital
part of the aircraft, fuse functioning and kill subject to being hit and fuse functioning. Pjh is hit on
the j-th vital part, in case of DA fused ammunition and landing in the vicinity zone, in case of
VT-fused ammunition. DA fuse ammunition is direct attack ammunition which hits the targets and
explodes. Variable time fuse ammunition (VT fuse ) explodes in close vicinity of the target. Vicinity
zone is a region around the aircraft in which once shell lands, will explode with some probability.
Probability of functioning is the design parameter of the fuse and varies from fuse to fuse.
Thus probability of kill of the aircraft, subject to kill of its vital parts, is given as
d id id
Pk /hf = 1 – 1 – Pk1,1 1 – Pk2,1 1 – Pk3,1 i (1– P4
P4
k,1 k,2 ) d1 – Pk5,1 Pk5,2 i ...(6.2)
where Pkj,i is the kill probability of the i-th redundant part of j-th vital part. In this relation Pkj,i is
same as Pk j in relation (6.1), except that subscript i has been added and bar over P has been removed.
This relation has been explained in chapter three in details.
If A and E are respectively the angles in azimuth and elevation of the aircraft, measured from
the ground based radar, then the direction cosines of line GO (line joining gun position and centre
of the aircraft) are given by
Let co-ordinates of point Pi, which are apexes of the triangular elements be (ui, vi, wi), in terms
of O-UVW co-ordinate system. In order to achieve our goal (to find projection of aircraft on N-plane)
we first find projection of a typical point P, of the aircraft on the N-plane. It is to be noted that point
P is a function of time t. Thus co-ordinates of a typical point P at the aircraft in terms of fixed
co-ordinate system (x, y, z) are obtained by linear transformations,
xp = x0 + l 1. u p + l 2 . v p + l3 . wp
yp = y0 + m1. up + m 2 . v p + m3 . wp ...(6.3)
zp = z0 + n1. up + n 2 . v p + n3 . wp
where subscript p denotes co-ordinates of point P and (x0, y0, z0) are the co-ordinates of the
centre O of the aircraft. The direction cosines of line GP (point P is different from the point O, which
is the centre of the aircraft), are
l p = xp /GP
mp = yp /GP ...(6.4)
np = zp /GP
and the angle between θ the lines GP and GO is
θ = cos (l0 l p + m0 m p + n0 n p )
–1
...(6.4a)
Simulation Model for Aircraft Vulnerability 135
(
1– n02 ) (1– n02 )
0
The point Q at which the line GP (produced, if necessary) meets the N-plane satisfies the equations,
GQ = GO/cosθ ; OQ = GO tanθ
Since line GQ is nothing but the extension of line GP, therefore co-ordinates of the point Q in
the GXYZ frame turns out to be
xq = GQ . l p , yq = GQ . m p , zq = GQ . n p
where (lp , mp , np) have been derived in equation (6.4).
N-plane
T-axis(+)
Z-axis
O Q
P
S-axis(+)
G
Y(–) X-axis
Thus, the direction cosines of the line OQ (which is line on the N-plane) with respect to the
G-XYZ frame are
d
l q = x q – x 0 /OQ i
mq = ( yq – y0 ) /OQ
nq = (z q )
– z0 /OQ
tq = OQ . sin φ
where φ is the angle, line OQ makes with OS axis in N-plane and is given by,
( )
φ = cos lq . ls′ + mq . ms′ + nq . ns′ n
where (ls′ , ms′ , ns′ ) are the direction cosines of the OS-axis with respect to G-XYZ frame. Derivation
of (ls′ , ms′ , ns′ ) is given in appendix 6.1.
1 1 s 2 t 2
Ph = ∫∫Fq 2 σ 2 σ 2 dsdt
exp
2πσ s σt
− + ...(6.6)
s t
where Fq is the projected region of the component triangles over N-plane, σs, σt, are the standard
deviations along OS and OT axis computed from system errors of the gun in the azimuth and elevation
plans respectively.
the AD guns (fitted with DA fuse ) to a given aircraft, following criteria have been adopted in the model.
(a) Considering the actual terminal velocity, mass and calibre of shell, penetration in the vital
as well as non-vital parts has been calculated.
(b) If the energy released by a shell is greater than the energy required (Table 6.1) to kill a
vital part, then probability of kill of that part is taken as 1.0 and probability of kill of aircraft
is calculated as per the Table 6.2.
DA fuse has an in-built delay mechanism. This is due to the fact that , shell first has to penetrate
the target and then explode. Kinetic energy of the projectile is responsible for the penetration of the
projectile in the vulnerable parts. When a projectile penetrates a surface, its velocity and mass both
reduce and thus total kinetic energy reduces. This is due to the resistance offered by the target. In
the case of aircraft (and same is true in case of other targets too) damage to a vital part is caused
due to the remaining energy which is left with the projectile after crossing the outer skin. Remaining
velocity of projectile after penetration in the vulnerable part is given by [75]
( )
Vr = V – xρD0 R 2V sin α / m cos θ ...(6.7)
where
V = normal striking velocity,
Vr = remaining velocity,
ρ = density of target,
α = nose cone angle,
θ = striking angle of projectile,
D0 = thickness of the target,
R = radius of projectile,
m = mass of projectile.
In the present model it has been assumed that if shell hits at any of the vital or non-vital
(rest of aircraft body) part, its kill probability depends on the factor, whether shell penetrates the aircraft
body or not. The kill probabilities given in Table 6.2 are for the K-kill (Chapter three) of aircraft,
subject to kill of vital part, when a small calibre impact fused high explosive projectile (in the present
case 23 mm shell) hits it. If shell is of higher calibre which release energy Q then these kill probabilities
are multiplied by factor
Pkj/ i = Pk / h Fk / i . E f
j j
...(6.8a)
where Pkj/ h is the probability of kill given by equation (6.17), and factor Fkj/ i is
138 System Modeling and Simulation
given in the Table 6.2. It is to be noted that factor Ef is multiplied with the total kill probability only in the
case of DA fused ammunition. In equation (6.8a), j is not the dummy index but indicates j-th vital part.
1.0
(4.5, 0.8)
Probability of fused functioning
0.8
0.6
P (r ) f
0.4
0.2
(6.5, 0.0)
0
2 4 6 8
Miss distance in meters (r)
7
where C = ∫ Pf (r ) dr
0
and
Pf (r) = 0.8, for r ≤ 4.5
= 0.4r + 2.6, for 4.5 ≤ r ≤ 6.0
= –0.2r + 1.4, for 6.0 ≤ r ≤ 6.5
= 0.0, for r > 6.5
where r being the normal distance from the surface of the target.
Probability of fuse functioning in the above equation is for a typical fuse (Fig. 6.3) and may be
different for different types of fuses. In case of proximity fused ammunition shell can land anywhere
in the vicinity region of the aircraft.
Vicinity region of a target is a region around it, so that if a shell lands within it, its fuse functions.
Probability of fuse functioning depends on the normal distance from the surface of the target. Let ms
be the maximum distance at which there is a probability that fuse will function. This distance is called
miss distance i.e., the distance beyond which if shell lands, fuse will not function. In fact, vicinity region
has an irregular shape depending on the reflected signals coming from different parts of the aircraft.
To model such a region mathematically is a tedious task. Thus we have assumed that this region is
of cylindrical shape with its axis, same as the axis of the aircraft. Consider a cylinder with axis along
the axis O-U of the aircraft, whose radius is Rv = Ra + ms and length 2(l0 + ms) mid-point of cylinder,
being the geometric centre of aircraft, Ra and ms being the radius of the aircraft and miss-distance.
Probability of landing and fuse functioning of projectile in terms of fixed co-ordinate system,
around aircraft is
1
PL f = Pf ( r – Ra ).exp{– ( z ) 2 + ( y ) 2 + ( z ) 2 / 2σ 2 } d x d y d z ...(6.10)
( )
3
2πσ
where
x = x − x0 ,
y = y − y0 ,
z = z − z0
and Pf (r – Ra) is the probability of fuse functioning. Ra is value of r at the surface of aircraft, r
being the radial distance from the u-axis. Relation (6.10) is a simple extension of two dimensional
Gaussian distribution to three dimensional case. Converting ( x , y , z ) co-ordinate, to moving
co-ordinates (u, v, w) with the help of linear transformations (inverse of equations 6.3) one gets,
1
. (u 2 + v 2 + w2 ) FG x , y, z IJ
PLf = P
( 2π .σ)3 f
(r – R a) exp −
2σ 2
J1
H u, v, w K du . dv . dw
140 System Modeling and Simulation
where
l1 l2 l3
FG x , y, z IJ m1 m2 m3
J1
H u, v, w K =
n1 n2 n3
is the Jacobian, used for transformation of co-ordinates from one set of co-ordinate system to other
set of co-ordinate system. Transforming above equation to cylindrical co-ordinates (r, θ, ζ) one gets,
1 r 2 + ζ2
P Lf = .
Pf (r – Ra ) exp –1/ 2 J .r . dζ.dr . dθ ...(6.11)
( ) σ2
3
2πσ2
The probability of kill of j-th vital part due to shell landing at a typical point P (r,θ,ζ )of vicinity shell is,
where Pk j/ h is the kill probability of j-th part subject to a hit on it (6.8a). Since the shell can land
any where in the vicinity region, the cumulative kill probability Pkj/ i of j-th part, due to i-th (say)
shell landing anywhere in the vicinity of the aircraft is obtained by integrating (6.12) over the whole
vicinity shell i.e.,
RL Rv 2π ( l0 + ms )
Expression for Pkj/ h will be derived in coming sections. Evaluation of Pkj/ h depends on the kill
criteria. In integral (6.13), limit RL is a typical distance from the aircraft such that if shell explodes
between Rν and RL, damage is due to blast as well as fragments, otherwise it is only due to fragments.
In the next section, evaluation of parameter RL has been discussed.
6.4.1 Determination of RL
The estimation of RL can be done on the basis of critical ‘impulse failure criterion’(page 134 of [14]).
This criterion essentially states that structural failure under transient loading can be correlated to critical
impulse applied for a critical time duration where the latter is assumed to be one quarter of the natural
period of free vibration of the structure. This critical impulse can be expressed as:
I 0 = ( ρ /E )1/ 2 . t . σ y ...(6.14)
where
E = Young’s modulus of the target material,
ρ = density of material,
t = thickness of the skin,
σy = dynamic yield strength.
Simulation Model for Aircraft Vulnerability 141
td 980 1 + ( z / 0.54 ) 10
=
W 1/ 3 1 + ( z / 0.02)3 1 + ( z / 0.74 ) 6 1 + (z / 6.9 ) 2
I = ∫ Pr dt
0
142 System Modeling and Simulation
where t d is the duration of the shock pulse. If we assume the reflected pressure pulse (which is
quite fair approximation) to be a triangular pressure pulse, then
Pr . td
I =
2
Taking the dimensions of the panel as a and b, the natural frequency (ω) of fundamental mode is [66]
Et 3
ω = π (1/a + 1/b )
2 2 2
12(1 – v 2 ) ρ
where
E = Young’s modulus,
ρ = density,
t = thickness,
ν = Poisson’s ratio.
Then the natural time period T, of panel is T = 2π/ω.
Keeping in view the above relations, we can simulate the value of z for which I > IC. This simulated
value of z will be equal to RL.
It has been assumed that an aircraft is assumed to be killed with some probability, if at least one
of its vital parts is killed. A vital component is treated as killed if the remaining energy of the fragment
penetrating the vital part is more than the energy required to kill it (KILL CRITERION). It is quite
possible that a fragment may not have sufficient energy required to kill a vital part. Singh and Singh
[66] have assumed that if at least k number of fragments whose total kinetic energy is equal to the
required energy Er (Table 6.1), should penetrate a part in order to achieve its kill. Thus probability
i
of kill Pkm , that at least k number of fragments, whose mass is greater than m, should penetrate is
given by Poisson’s law as,
k −1
(mr ) N − mr
i
Pkm = 1− ∑ ∑ e ...(6.16)
N =0 N!
where mr is the average number of fragments penetrating the component. In this assumption
there is a weakness. If a fragment hits but does not penetrate, may be due to the fact that
it does not have sufficient energy for doing so, then in such a situation it will not cause any
damage to the component. This mean fragment should have sufficient energy to damage the
component. If a fragment does not have sufficient energy required for penetration in aircraft’s
skin, it will be reflected from its surface. Even if there are ten such fragments will not be
able to damage the aircraft. Keeping this in mind this criterion later was modified by another
criterion, known as energy criterion and is being discussed below.
conical zones, so that five degree is the apex angle of each zone. These zones are called uniform
because average number of fragments per unit solid angle is constant with in a particular zone.
Fig. 6.5: Pictorial view of a shell exploding in the near vicinity of aircraft.
(Straight lines emerging from shell centre are fragment zones.)
Since the size of the shell is very small as compared to that of aircraft, it can be assumed that
fragments are being ejected as if they are coming from a point, which is centre of the shell (Fig. 6.5).
k
We define, f i , i +1 as the region, which is the intersection of two solid cones, with vertex at point Pk and
whose slant surfaces make angles αi, αi+1 respectively with the axis of the shell, positive direction being
k
towards the nose of the shell. Let ni , i +1 be the total number of fragments of mass greater than m, in the
angular zone n ik,i +1 . Fragments per unit solid angle in the z ik,i +1 -th angular zone can be given as
n ki, i +1
f ik, i +1 = ...(6.18)
2π(cos αi – cos αi +1 )
It is to be noted that angles αi, αi+1 are the angles of fragment cones in static mode i.e., when
shell has no velocity. Distribution and number of fragments after explosion are obtained experimentally
by exploding a static shell in air and measuring the distribution of fragments by collecting them in
cardboard cabin made for the purpose. But it is known that in actual practice shell is moving with
reference to fixed co-ordinate system and these angles have to be evaluated in that mode. Let these
angles in dynamic mode be α' j . α' j +1 (see relation 6.22).
k
Let ω ik,i+1 be the solid angle subtended by the component in z i, i +1 the angular zone. The total
number of fragment hits on a component from this zone will be given by
nz
Nk = ∑ω i =1
k
i , i +1
. f ik,i +1 ...(6.19)
Simulation Model for Aircraft Vulnerability 145
These are the average number of fragments hitting a vital part, which is used in equation (6.16)
and is denoted by a symbol mr.
The solid angle subtended in the angular zone z i , i + 1 (here we have omitted the superscript k
as these equations hold good for all k’s) by a component at the centre of gravity of the shell is
determined by the intersecting surface of the component and the angular zone z i , i + 1 , and can be
given by (Fig. 6.6)
ωi.i+1 = ∑ δw ...(6.20)
Ai , i +1
where
cos θ δA
δw = ...(6.20a)
RA2
where
Ai,i +1 = intersecting surface of the component and the zone zi,i + 1 which will differ in
stationery and dynamic cases.
A = a small elemental area on the vital component whose surface area is Ai,i+1.
R A = distance between centre of gravity of the shell and mid point of δA.
θ = angle between RA and normal to the surface at the mid point of δA.
Value of δw is evaluated in equation (6.20a). Following is the example to evaluate the solid angle
subtended by a component of the aircraft in the different angular zones of the PF-shell, when the
shell bursts at any arbitrary point Cs in the vicinity region of the aircraft. Using this method kill for
any component of the aircraft can be obtained, having well defined surface.
Let the PF-fused shell bursts at a point Cs (also it is centre of gravity of the shell) in the vicinity
region of the aircraft say at time t = 0. Let the co-ordinates of the centre of gravity of the shell at
time t = 0 is ( x s , y s , z s ) and velocity is vs in the direction of (ls, ms, ns) which is the direction of
its axis with respect to I-frame, fixed in space.
Further let at the time of burst, ( x a , y a , z a ) are the co-ordinates of the centre of the aircraft
which is also the origin of the II-frame and let,(li, mi , ni), i = 1, 2, 3 are the direction cosines
of the aircraft’s axes (i.e., axes of the second frame) with respect to frame-I and this aircraft
(or frame-II) is moving with velocity Va in the direction (lv, mv, nv) in the frame-I.
If co-ordinates of centre of shell at the time of burst (t = 0) with respect to frame-II of reference
are (u s , v s , w s ) , then transformation from one system of co-ordinates to other is given as
xs = us . l1 + vs . l2 + ws . l3 + xa
ys = us . m1 + vs . m2 + ws . m3 + ya
zs = us . n1 + vs . n2 + ws . n3 + za
146 System Modeling and Simulation
us = xs . l1 + ys . m1 + zs . n1 + xa
v s = xs . l2 + ys . m2 + zs . n2 + ya
ws = xs . l3 + ys . m3 + zs. n3 + za (6.21)
Let us assume that PF-fused shell bursts in stationary position with reference to frame-I αi, αi + 1
are the angles which the boundaries of the conical angular zone of fragments zi , i +1 make with the
positive direction of the shell axis and V f i , V f ( i + 1) are the corresponding velocities of the fragments
of these boundaries.
When shell bursts in dynamic mode, the directions and velocities of fragments, as
observed in frame-I will be,
α′i = tan –1 (V2 / V1 )
d
Vfi¢ = V 12 + V 22 i
1/ 2
...(6.22)
where
V1 = Vs + V f i cos ( α i )
V2 = V fi .sin ( αi )
Fragments emerging from the point Cs in an angular zone zi, i + 1 will be confined in a cone making
angles α′i and α′i + 1 respectively with the axis of the shell. Intersection of this cone with component
is say an area P1, P2, P3, P4. Divide the surface enveloping P1, P2, P3, P4 into a finite number of
rectangular areas δA = δl . δb (say) (Fig. 6.6) where δ l and δb are dimensions of the rectangular
element.
d i
If a typical point P whose co-ordinate with respect to frame-II are u p , v p , w p is the middle point
of area δA, then solid angle of area δA subtended at Cs and angular zone to which it belongs is
determined by simulation as follows.
Co-ordinates of point P with reference to the fixed frame, at time t after burst are
xpt = xp + Va lv .t
xpt = yp + Va mv .t ...(6.23)
zpt = zp + Va nv .t
If θ is the angle between shell axis and the line Cs – Pt where point Pt ( x pt , y pt , z pt ) being the position
of P at time t, then first step is to determine the angular zone α′i, α′i+1 in which θ lies.
Fragment may come to the point Pt from angular zone zi, i + 1 with velocity F ′f i F f′ (i + 1) depending
upon is θ close to.
Distance travelled by the fragment along the line Cs – Pt in time t
Df = V f′ . t ...(6.24)
d i
where Vf¢ = selected V ft V f ( i +1) . This means Vf¢ takes the value Vf¢ , or Vf¢ (i +1)
¢ ¢
′ ′
subject to the condition that it is near to αi or αi +1 . Actual distance between Cs and Pt is
2 1/ 2
(
Ds = Σ xs – x pt
) ...(6.25)
From equations (6.24) and (6.25) we simulate time t, such that Df = Ds for confirmed impact.
The velocity of impact to the aircraft, Vstrike (relative impact velocity) is,
( )
2 1/ 2
V strike = V ' f + Va2 – V f′ Va cos β ...(6.26)
148 System Modeling and Simulation
where β is the angle between the positive direction of aircraft’s velocity vector and fragment’s velocity
vector. Penetration for strike velocity Vs can be obtained from experimental data [14]. Penetration
depends on angle of impact. If θ is the angle of impact, then
V θ = V0/cosθ ...(6.27)
where Vθ is the required velocity of impact at angle θ when velocity at zero angle of impact is given.
If velocity of impact Vstrike is such that it is more than some critical velocity then the solid angle
δω, subtended by the small rectangular element, in the angular zone z i , i +1 , given by
dA .|cos θ |
δω = …(6.28)
Ds2
is added to the relation (6.20).
a4
( ) 1
a2
mr = ms –10 a1 kD 3 ms2 msa3 a5
Vs
cos θ
b4
( ) 1
b2
Vr = Vs –10b1 kD 3 ms2 msb3 Vsb5 ...(6.29)
θ
cos
where values of ai and bi for duralumin is given in Table 6.3. ms is striking mass in grains, Vs is
striking velocity in ft/sec, D thickness in inches, k is the shape factor for the projectile.
Shape factor k = .0077 spherical fragments
= .0093 cubical fragments
Pkj = 1– ∏ 1– Pk / i
j
...(6.30)
i =1
Further the aircraft can be treated as killed if at least one of its vital parts is killed. Thus in this
case the cumulative kill probability for the aircraft as a whole can be given as:
y
Pk = 1– ∏ 1– Pk
j
...(6.31)
i =1
where factor Pk j/ l is the probability of kill of aircraft subject to a vital component kill (Table 6.2 and
eqn. (6.8a)). In case the j-th part of the aircraft has yi redundant parts and P is given by
y y
Pk = 1 − ∏ 1– ∏ 1– Pk j/ r ...(6.31a)
i =1
j =1
where Pk j/ r denotes the kill probability of the r-th redundant part of the j-th vital part.
Each vital part is covered by some of the above mentioned triangular faces of the structure through
which it can get lethal hits. Such triangles have been used to estimate kill probability of that vital
component. Here pilot 1 and pilot 2 are nonredundant but fuel tanks and engines has been taken as
redundant parts.
– Data of air defence gun
An Air Defence twin barrel gun with DA/VT fused ammunition is considered. Data of a typical
air defence gun is given in Table 6.5 and critical and uncritical energy required to kill vital parts of
the present aircraft has been given in Table 6.4.
Table 6.5: Gun data
Parameters Value
System error 3 m rad
Number of barrels 2
Firing rate 5 rounds/gun barrel
Probability of DA fuze functioning 0.99
Probability of proximity fuze functioning see Fig.7.3
Time of continuous firing of gun 3 seconds
Maximum range of gun 5000m
Minimum range of gun 500m
Maximum detecting range 10000m
Results of the model are computed and are shown in Table 6.6 for DA fused ammunition and
in Table 6.7 for PF fused ammunition. First column gives cumulative kill probability of aircraft and
other columns give kill probabilities of individual components. Similarly in Fig. 6.9, variation of kill
probability of an aircraft vs number of rounds fired is shown when aircraft is engaged at range two
kilo meters.
Table 6.6: Variation of probability of kill of aircraft and its five vital parts
due to DA-fused ammunition (Engagement range = 2000m)
Rounds Aircraft Avionics Pilot 1 Pilot 2 Fuselage Wing fuel Engine1 Engine 2
fuel tank tank
2 0.0122 0.0008 .0100 0.0015 0.0009 0.0152 0.0010 0.00100
Contd...
Simulation Model for Aircraft Vulnerability 151
It is seen that kill due to PF fused is much higher than that DA fused ammunition shell and increased
with the number of rounds. Figure 6.10 gives variation of kill vs. Range. It is seen that kill decreases
with higher opening range, for thirty rounds.
Figure 6.8 gives computer output of an aircraft generated by combining the triangular data from the
drawings.
Table 6.7: Variation of probability of kill of the aircraft and its five vital parts
due to VT-fused ammunition at range 2.0 km
Rounds Aircraft Avionics Pilot 1 Pilot 2 Fuselage Wing fuel Engine 1 Engine 2
fuel tank tank
2 0.0293 0.0193 0.0044 0.0037 0.0384 0.0441 0.0195 0.0196
4 0.0605 0.0380 0.0085 0.0071 0.0760 0.0870 0.0388 0.0389
6 0.0947 0.0564 0.0126 0.0105 0.1131 0.1290 0.0583 0.0585
8 0.1312 0.0747 0.0167 0.0141 0.1498 0.1701 0.0780 0.0783
10 0.1697 0.0928 0.0210 0.0177 0.1858 0.2103 0.0980 0.0983
12 0.2097 0.1106 0.0253 0.0214 0.2213 0.2495 0.1181 0.1185
14 0.2509 0.1281 0.029 0.0252 0.2562 0.2877 0.1384 0.1389
16 0.2928 0.1454 0.0343 0.0292 0.2904 0.3249 0.1588 0.1595
18 0.3350 0.1624 0.0389 0.0332 0.3239 0.3610 0.1795 0.1802
20 0.3772 0.1791 0.0436 0.0373 0.3567 0.3961 0.2003 0.2011
22 0.4191 0.1956 0.0484 0.0416 0.3889 0.4302 0.2212 0.2221
24 0.4604 0.2117 0.0533 0.0460 0.4203 0.4632 0.2424 0.2434
26 0.5010 0.2275 0.0583 0.0505 0.4510 0.4951 0.2637 0.2648
28 0.5404 0.2430 0.0634 0.0552 0.4809 0.5260 0.2851 0.2863
30 0.5785 0.2581 0.0686 0.0600 0.5100 0.5557 0.3067 0.3080
152 System Modeling and Simulation
0.6
DA
VT
0.4
CKP
0.2
0
2 8 14 20 26
No. of rounds
Fig. 6.9: Variation of CKP vs. number of rounds (Engagement range 2000m).
0.8 DA
0.6 VT
CKP 0.4
0.2
0
2000 2500 3000 3500
Range
APPENDIX 6.1
Let (x0, y0, z0 ) be the centre of the N-plane forming a right handed system of axis, S-axis,
T-axis and OG axis as shown in Fig. 6.2. The line GO is perpendicular to the ST plane with direction
cosines (–l0, –m0, –n0) where
x0
l0 = cosA .cosE =
| GO |
y0
m0 = sin A . cosE = ...(B.1)
GO
x0
n0 = sin E =
GO
S-axis which will lie in the so-called azimuthal plane will be normal to the elevation plane i.e.,
normal to the plane GOO′ where O′ is the projection of point O on GXY plane.
Equation of the plane GOO′ is
m −l0
0
, ,0 ...(B.6)
2
1 − n02
1 − n0
154 System Modeling and Simulation
lt x0 + mt y0 + nt z0 = 0 ...(B.7)
Solving equations (B.7), (B.8) and (B.9) for lt, mt, nt we get the directions cosines of OT axis
as
−n l mn
0 0
, − 0 0 , 1 − n02 . ...(B.10)
2
1 − n0 1 − n02
Simulation Model for Aircraft Vulnerability 155
APPENDIX 6.2
l = (u c − u g ) D s
m = ( vc − v g ) D s
n = ( wc − wg ) Ds
156 System Modeling and Simulation
Ds = ∑ (uc − u g ) 2
DCs’ (a, b, c) of the normal to a triangular surface co-ordinates of whose corners are (ui, vi, wi),
i = 1, 2, 3 can be obtained as
v1 w1 1
A = v2 w2 1
v3 w3 1
u1 w1 1
B = u2 w2 1
u3 w3 1
u1 v1 1
C = u2 v2 1
u3 v3 1
A
a=
A + B2 + C 2
2
B
b=
A + B2 + C 2
2
C
c =
A + B2 + C 2
2
Let α is the angle between line GC and normal to the triangular face then
α = cos ( al + bm + cn )
−1
Let (sik, tik), k = 1, 3 are the co-ordinates of the corners of projection of i-th triangle of the aircraft
over N-plane with respect to point G and (sjk, tjk), k = 1, 3 are the co-ordinates of the projection
of j-th triangle over N-plane. These triangles are further subdivided into smaller rectangular meshes,
in order to assess the overlapped area.
Let ( s0 , t0 ) be the central point of a typical rectangle of the i-th triangle. If this point falls on
in the j-th triangle formed by the vertices (sjk, tjk), k = 1, 3 then it implies that this rectangle is overlapped
by the j-th triangular face. If it is not covered by j-th triangle, other triangles are tested. Similar test
is applied to all the rectangular elements of the i-th triangle. Thus if all the rectangular elements are
not covered by any other triangle, it implies that this triangle is not being overlapped by any of the
triangles and can be considered to find solid angle or hit probability. Let a rectangle with centre
c p ( s0 , t0 ) is overlapped by j-th triangle. In that case we have to check whether the j-th triangular
face is near to the source point of the projectile or the i-th triangular face is nearer. Whichever triangle
is nearer will be overlapping other. It can be done in the following steps.
(i) Let us define a III-co-ordinate system, OST, with origin at 0 and s-t plane being normal
to line joining projectile and centre of the aircraft.
(ii) Let direction cosines of the line, joining points c p ( s 0 , t 0 ) and source point of the projectile
G with respect to frame-III are (l3′ , m3′ , n′ 3 )
(iii) ( )
Let l1′ , m1′ , n1′ be the direction cosines of GC p with reference to frame-I. Thus
where (ls , m s , ns ); (lt , mt , nt ) and (lr , mr , nr ) are the direction cosines of the axes of
frame-III with respect to frame-I.
Thus and the DC’s of the line GC p with respect to frame-II say (l′2, m′2, n′2)
l′ 2 = l ¢1 l1 + m1¢ m1 + n1¢ n1
m ¢2 = l ¢1 l2 + m ¢1 m2 + n1¢ n2
n 2¢ = l ¢1 l3 + m¢1 m3 + n1¢ n3
(iv) Finally find the equation of the line with respect to frame-II as the line passes though some
point whose co-ordinate with respect to frame-II are known.
(v) Let the line meet the i-th triangular face at point O i with co-ordinates oi (uoi , voi , woi )
uoi = u g + l2′ R0
uoi = vg + m ¢2 R0
uoi = wg + n2¢ R0
158 System Modeling and Simulation
au g + bvg + cwg + di
R0 = –
al2′ + bm2′ + cn2′
12345678
12345678
12345678
12345678
12345678
7
12345678
12345678
SIMULATION OF QUEUING
SYSTEMS
1. A.K. Erlang was the first person to study the problem of telephone networks. By studying a village telephone
exchange he worked out a formula, now known as Erlang’s formula, to calculate the fraction of callers attempting to
call someone outside the village that must wait because all of the lines are in use. Although Erlang’s model is a simple
one, the mathematics underlying today’s complex telephone networks is still based on his work.
He was born at Lønborg, in Jutland, Denmark. His father, Hans Nielsen Erlang, was the village schoolmaster and
parish clerk. His mother was Magdalene Krarup from an ecclesiastical family and had a well known Danish
mathematician, Thomas Fincke, amongst her ancestors. He had a brother, Frederik, who was two years older and
two younger sisters, Marie and Ingeborg. Agner spent his early school days with them at his father’s schoolhouse.
Evenings were often spent reading a book with Frederik, who would read it in the conventional way and Agner would
sit on the opposite side and read it upside down. At this time one of his favourite subjects was astronomy and he
liked to write poems on astronomical subjects. When he had finished his elementary education at the school he was
given further private tuition and succeeded in passing the Praeliminaereksamen (an examination held at the University
of Copenhagen) with distinction. He was then only 14 years old and had to be given special entrance permission.
Agner returned home where he remained for two years, teaching at his father’s school for two years and continuing
with his studies. He also learnt French and Latin during this period. By the time he was 16 his father wanted him to
go to university but money was scarce. A distant family relation provided free accommodation for him while he prepared
160 System Modeling and Simulation
made to model science of queues. The basic concept of queuing theory is the optimization of wait
time, queue length, and the service available to those standing in a queue. Cost is one of the important
factors in the queuing problem. Waiting in queues incur cost, whether human are waiting for services
or machines waiting in a machine shop. On the other hand if service counter is waiting for customers
that also involves cost. In order to reduce queue length, extra service centers are to be provided
but for extra service centers, cost of service becomes higher. On the other hand excessive wait
time in queues is a loss of customer time and hence loss of customer to the service station. Ideal
condition in any service center is that there should not be any queue. But on the other hand service
counter should also be not idle for long time. Optimization of queue length and wait time is the
object theory of queuing.
Let us see how this situation is modeled. First step is to know the arrival time and arrival
pattern of customer. Here customer means an entity waiting in the queue. One must know from
the past history, the time between the successive arrival of customers or in the case of machine
shop, the job scheduling. Also arrival of number of customers vary from day to day. On Saturdays,
number of customers may be more than that on other days. What is the probability that a customer
will arrive in a given span of time, is important to know.
In order to maximize the profit, the major problem faced by any management responsible for
a system is, how to balance the cost associated with the waiting, against the cost associated with
prevention of waiting. An analysis of queuing system will provide answers to all these questions.
However, before looking at how queuing problem is to be solved, the general framework of a queuing
system should be understood.
A queuing system involves customers arriving at a constant or variable time rate for service at
a service station. Customers can be students waiting for registration in college, aeroplane queuing
for landing at airfield, or jobs waiting in machines shop. If the customer after arriving, can enter
the service center, good, otherwise they have to wait for the service and form a queue. They remain
in queue till they are provided the service. Sometimes queue being too long, they will leave the
queue and go, resulting a loss of customer. Customers are to be serviced at a constant or variable
rate before they leave the service station. A typical queuing system is shown in Fig. 7.1.
for his University entrance examinations at the Frederiksborg Grammar School. He won a scholarship to the University
of Copenhagen and completed his studies there in 1901 as an M. A. with mathematics as the main subject and
Astronomy, Physics and Chemistry as secondary subjects.
Over the next 7 years he taught in various schools. Even though his natural inclination was towards scientific
research, he proved to have excellent teaching qualities. He was not highly sociable, he preferred to be an
observer, and had a concise style of speech. His friends nicknamed him “The Private Person”. He used his
summer holidays to travel abroad to France, Sweden, Germany and Great Britain, visiting art galleries and
libraries. While teaching, he kept up his studies in mathematics and natural sciences. He was a member of the
Danish Mathematicians’ Association through which he made contact with other mathematicians including members
of the Copenhagen Telephone Company. He went to work for this company in 1908 as scientific collaborator and
later as head of its laboratory.
Erlang at once started to work on applying the theory of probabilities to problems of telephone traffic and in 1909
published his first work on it “The Theory of Probabilities and Telephone Conversations”1 proving that telephone calls
distributed at random follow Poisson’s law of distribution. At the beginning he had no laboratory staff to help him, so he
had to carry out all the measurements of stray currents. He was often to be seen in the streets of Copenhagen,
accompanied by a workman carrying a ladder, which was used to climb down into manholes.
Simulation of Queuing Systems 161
Input Output
(Customers) (Customers)
Waiting line Service facility
Because of the growing interest in his work several of his papers were translated into English, French and German. He
wrote up his work in a very brief style, sometimes omitting the proofs, which made the work difficult for non-specialists
in this field to understand. It is known that a researcher from the Bell Telephone Laboratories in the USA learnt Danish in
order to be able to read Erlang’s papers in the original language.
His work on the theory of telephone traffic won him international recognition. His formula for the probability of loss was
accepted by the British Post Office as the basis for calculating circuit facilities. He was an associate of the British
Institution of Electrical Engineers.
Erlang devoted all his time and energy to his work and studies. He never married and often worked late into the night. He
collected a large library of books mainly on mathematics, astronomy and physics, but he was also interested in history,
philosophy and poetry. Friends found him to be a good and generous source of information on many topics. He was
known to be a charitable man, needy people often came to him at the laboratory for help, which he would usually give
them in an unobtrusive way. Erlang worked for the Copenhagen Telephone Company for almost 20 years, and never
having had time off for illness, went into hospital for an abdominal operation in January 1929. He died some days later on
Sunday, 3rd February 1929.
Interest in his work continued after his death and by 1944 “Erlang” was used in Scandinavian countries to denote the unit
of telephone traffic. International recognition followed at the end of World War.
162 System Modeling and Simulation
1. Calling source. or the population from which customers are drawn. Calling source may
be finite or infinite. When queue is so long that arrival of one more customer does not
effect the queue length, we call it infinite source of customers. A reverse of this situation,
when queue length is not long and incoming or outgoing of one-customer affects the queue;
we call it a finite source of customers.
2. The input or arrival process. This includes the distribution of number of arrivals per unit
of time, the number of queues that are permitted to be formed, the maximum queue length,
and the maximum number of customers desiring service.
3. The service process. This includes time allotted to serve a customer, number of servers
and arrangement of servers. Simplest case is single queue and single server.
e −λ λ x x = 0, 1, 2,...
f ( x) = Pr( X = x) = , ...(7.1)
x! λ > 0
E( X ) = λ
where λ is the average number of arrivals per unit time (1/τ), and x is the number of customers per
unit time. This pattern of arrival is called Poisson’s arrival pattern.
It is interesting to know that second assumption leads us to the result that inter arrival time T
follows an exponential distribution1, with the same mean parameter λ. To prove this, let us assume
T = time between consecutive arrivals. If an arrival has already occurred at time t = 0, the time to
the next arrival is less than t, if and only if there are more than one arrival in time interval [0, t].
This probability G(t) that inter-arrival time is less than t can be defined as,
∞
G(t ) = Pr(T < t ) = ∑ e − λ t (λt ) x / x !
x =1
Therefore
∞
∑e
x =1
− λt
( λt ) x / x! = 1 – e–λt
−λt
G(t ) = Pr(T < t ) = 1 − e
1. An Alternative proof
∆t
The probability of arrival of a customer during a very small time interval ∆t is . Hence probability of a customer
τ
∆t
not arriving during time ∆t is (1 – ). Now if
τ
h(t) = probability that the next customer does not arrive during the interval t given that the previous
customer arrived at time t = 0, and likewise.
h(t + ∆t ) = probability that the next customer does not arrive during the interval (t + ∆t ) given that the
previous customer arrived at time t = 0.
Since the arrival of customers in different periods are independent events (i.e., the queue has no memory), we
write
∆t
h(t + ∆ t ) = h(t ) . (1 − ) ...(7.1a)
τ
h (t + ∆ t ) − h (t ) h (t )
or = −
∆t τ
First equation of (7.1a) means, probability that customer does not arrive in interval (t + ∆ t) is equal to the
probability that he does not arrive in the interval t and also in the interval ∆ t.
Contd...
164 System Modeling and Simulation
d [G (t )] d (1 − e −λt )
g (t ) = = = λe −λt ...(7.2)
dt dt
Equation (7.2) is the exponential probability density function discussed earlier in chapter two.
Figure 7.2 gives plot of exponential density function and its cumulative distribution function.
The curve in Fig. 7.2(b) gives the probability that the next customer arrives by the time t, given
that the preceding customer arrived at time zero. It is appropriate to mention here that inverse of
inter arrival time τ is denoted by λ and is the average number of customers at the server per unit
time.
1 1
f (x) F(x)
0
t/ t/
(a) Exponential density function (b) Exponential distribution function
g (t ) =
τ
(
1 –t /τ
e ) ...(7.3a)
which is the distribution for inter arrival time.
Simulation of Queuing Systems 165
e −λ λ x e −1212 x x = 0, 1, 2,...
f (x) = Pr(X = x) = = ,
x! x! λ > 0
E ( X ) = 12 cars/hr
The distribution of the time between consecutive arrivals T is,
−12t
g (t ) = 12e ; t > 0
1
E (T ) = hr between arrivals.
12
Third assumption (steady state) means queuing system has been operating long enough to be
independent of the initial state of the system and is independent of time. That is, the system has reached
a state of equilibrium with respect to time. The distribution of the number of arrivals per unit time
and the distribution of the service time do not change with time.
q1 (t + ∆t ) = (probability that no arrival takes place between time zero and t) . (probability
that single arrival takes place during time ∆t) + (probability that single arrival
takes place between time zero and t) . (probability that no arrival takes place
during time ∆t).
∆t ∆t
= f (t ) . + q1 (t ) . (1 − )
τ τ
q1 (t + ∆ t ) − q1 (t ) 1
= [ f (t ) − q1 (t )]
∆t τ
dq1 1
= [ f (t ) − q1 (t )]
dt τ
166 System Modeling and Simulation
t −t / τ t
q1 (t ) = e = f (t )
τ τ
Now we extend the above logic for two customers in the queue.
q 2 (t + ∆ t ) = (probability that one arrival takes place between time zero and t).
(probability that single arrival takes place during time ∆t)
+ (probability that two arrivals take place between time zero and t).
(probability that no arrival takes place during time ∆t).
∆t ∆t
= q1 (t ) . + q2 (t ) 1 −
τ τ
q2 (t + ∆t ) – q2 (t ) 1
= [q1 (t ) – q2 (t ) ]
∆t τ
When limit ∆ t → 0, this equation becomes,
dq2 1
= [ q1 (t ) − q2 (t )]
dt τ
Above equation can be integrated as,
2
1t
q2 (t ) = f (t )
2! τ
We generalize this logic for k customers to arrive between time zero and t as,
k
1t
q k(t) = f (t ) ...(7.4)
k ! τ
where ρ is called the utilization factor of the service facility. This is also the average number of
customers in the service facility. Thus probability of finding service counter free is
(1 – ρ ) ...(7.6a)
That is there is zero customer in the service facility. Let Pn(t) be the probability of exactly n
customers being in the system at time t. Let ∆t > 0 be a small interval of time. The probability of
one customer arriving and no customer departing during the interval ∆t is
λ . ∆ t . (1 − µ ∆ t)
Similarly, probability of one customer arriving and one customer leaving during the interval ∆t is
(λ . ∆ t) . (µ . ∆ t)
The probability of no customer arriving and one customer leaving is
(1 − λ . ∆ t) . (µ . ∆t)
168 System Modeling and Simulation
+ Pn +1 (t ) . (1 − λ∆t ) . (µ∆t )
+ Pn −1 (t ) . (λ∆t ) . (1 − µ∆t )
From the above equation one gets,
Pn (t + ∆t ) − Pn (t )
= Pn (t ) . (λµ∆t ) + Pn (t )[−λ − µ + λ . µ∆ t)
∆t
+ Pn +1 (t ) . (1 − λ∆t ) . (µ )
+ Pn −1 (t ) . (λ ) . (1 − µ∆t )
= µ Pn +1 ( t ) + λ Pn −1 (t ) − Pn (t )( λ + µ )
+λ . µ∆ t [2 Pn (t ) − Pn+1 (t ) − Pn −1 (t )]
Taking limits of both sides of this equation when ∆t tends to zero,
dPn (t )
= µPn +1 (t ) + λPn −1 (t ) − (λ + µ ) . Pn (t ) ...(7.7)
dt
This equation holds for all n > 0. When n = 0, the contributions made by the term Pn – 1 would
be zero. Therefore
dP0 (t )
= µP1 (t ) − λ.P0 (t ) ...(7.8)
dt
If ρ < 1, after the passage of sufficiently long time the queue would reach an equilibrium, and
Pn(t) would converge to a constant. Thus its derivative can be put equal to zero in equilibrium condition.
The equations (7.7) and (7.8) become,
λ λ
Pn+1 (t ) + Pn −1 (t ) − + 1 . Pn (t ) = 0
µ µ
Pn +1 (t ) = (1 + ρ) . Pn (t ) − ρPn −1 (t ) , for n ≥ 1 ...(7.9)
λ.
P1 (t ) = P0 (t ) = ρP0 (t ) ...(7.10)
µ
Simulation of Queuing Systems 169
Now since ∑P
n=0
n = 1, therefore
∞
P0 ∑ ρ n = 1
n =0
1
or P0 = 1
1− ρ
since ρ < 1, therefore P0 = 1 − ρ
which means there is no body (zero person in the system) in the system (queue plus server) and
service counter is free, which is the same result as in equation (7.6a).
We define average number of customers at time t, in the system as LS . Here bar over L depicts
average length of queue and subscript S is for system. Similar terminology will be used for other
symbols in the following sections. LS is given by
∞
ρ
LS = ∑ nPn = P ∑ nρ
0
n
=
1− ρ
n= 0
λ ρ
Thus LS = = ...(7.12)
µ − λ (1 − ρ)
These and similar other statistics about the queue are called the operating characteristics of the
queuing system.
Average number of customers in the queue LQ is same as expected number in the system – the
expected number in the service facility i.e.,
λ λ λ2 ρ2
LQ = LS − ρ = − = = ...(7.16)
µ−λ µ µ(µ − λ ) (1 − ρ)
Average time a customer spends in the system is denoted by W S , and is equal to expected number
of customers in the system at time t, divided by number of customers arrived in unit time i.e.,
λ .1 1
WS = = ...(7.17)
µ − λ λ (µ − λ )
( )
Average time a customer spends in the queue WQ is same as average time a customer spends
in the system – average time a customer spends in the server i.e.,
1 λ
WQ = WS − =
µ µ (µ − λ )
= average time a customer spends in the queue. ...(7.18)
Similarly few other parameters can be defined as follows.
Now expected time T to serve a customer is given by,
1
E (T = t ) = =υ
µ
and time for which server remains idle in t seconds is given by (1– ρ) t/ υ , thus
Probability that the time in the system is greater than t is given by,
P(T > t ) = e −µ (1−ρ )t ...(7.19)
Similarly probability of more than k customers in the system is,
k +1
λ
P( n > k ) = ...(7.20)
µ
Below, we give few examples to illustrate these statistics.
Example 7.2: In a tool crib manned by a single assistant, operators arrive at the tool crib at the
rate of 10 per hour. Each operator needs 3 minutes on the average to be served. Find out the loss
of production due to waiting of an operator in a shift of 8 hours if the rate of production is
100 units per shift.
Solution: Arrival rate (λ) = 10 per hour
Service rate (µ) = 60/3 = 20 per hour
Average waiting time in the queue
λ 10 1
(W )Q =
µ (µ − λ )
=
20(20 − 10)
=
20
hour
Simulation of Queuing Systems 171
100 2
∴ loss of production due to waiting = × = 5 units
8 5
Example 7.3: At the ticket counter of football stadium, people come in queue and purchase tickets.
Arrival rate of customers is 1/min. It takes at the average 20 seconds to purchase the ticket.
(a) If a sport fan arrives 2 minutes before the game starts and if he takes exactly 1.5 minutes
to reach the correct seat after he purchases a ticket, can the sport fan expects to be seated
for the tip-off ?
(b) What is the probability the sport fan will be seated for the start of the game?
(c) How early must the fan arrive in order to be 99% sure of being on the seat for the start
of the game?
Solution: (a) A minute is used as unit of time. Since ticket is disbursed in 20 seconds, this means,
three customers enter the stadium per minute, that is service rate is 3 per minute.
Therefore,
λ = 1 arrival/min
µ = 3 arrivals/min
1
(W )
S = waiting time in the system =
µ–λ
=
1
2
The average time to get the ticket and the time to reach the correct seat is 2 minutes exactly,
so the sports fan can expect to be seated for the tip-off.
(b) This is equivalent to the probability the fan can obtain a ticket in less than or equal to half
minute.
P (T < 1/2) = 1 – P(T > 1/2)
λ
1 −µ (1 − ) t
P T < = 1 − e µ
2
1 1
–3(1 − ) . ( )
= 1− e 3 2
= 1 − e −1 = 0.63
(c) For this problem, we need to determine t such that, probability of arrival of fan in the
playground just before the start is 0.99, or probability that he does not arrive in time is 0.01.
λ
1 −µ (1 − ) t
P T < = 1 − e µ
2
λ
−µ (1 − )t
µ
But P(T > t ) = e
1
− 3(1 − )t
⇒ =e 3
= 0.01
⇒ −2t = ln (0.01)
∴ t = 2.3 minutes
172 System Modeling and Simulation
Thus, the fan can be 99% sure of spending less than 2.3 minutes obtaining a ticket (awaiting
for the purchase of ticket). Since it can take exactly 1.5 minutes to reach the correct seat after
purchasing the ticket, the fan must arrive 3.8 minutes (2.3 + 1.5 = 3.8) early to be 99% of seeing
the tip-off.
Example 7.4: Customers arrive in a bank according to a Poisson’s process with mean inter arrival
time of 10 minutes. Customers spend an average of 5 minutes on the single available counter, and
leave. Discuss (a) What is the probability that a customer will not have to wait at the counter, (b)
What is the expected number of customers in the bank, (c) How much time can a customer expect
to spend in the bank.
Solution: We will take an hour as the unit of time. Thus
λ = 6 customers/hour,
µ = 12 customers/hour.
(a) The customer will not have to wait if there are no customers in the bank. Thus
λ
P0 = 1 −
µ
= 1 − 6 /12 = 0.5
(b) Expected number of customers in the bank are given by
λ 6
LS = = =1
µ−λ 6
(c) Expected time to be spent in the bank is given by
1 1
WS = = = 1/6 hour = 10 minutes
µ−λ 12 − 6
So far, we have discussed a system consisting of single queue and single service counter. A
system may have a single queue and multiple service counters as shown in Fig. 7.4. So far, we assumes
that arrival at tern follows Poisson’s distribution and time for servicing follows an exponential
distribution. But these are not the only patterns of arrival. Arrival may be totally irregular, with time
as well as with day. For examples in banks, crowd is much more on Saturdays than on other days,
and thus Poisson’s pattern fails. Although mathematical models have been developed for different arrival
patterns, yet generally industrial and business units adopt simulation techniques for this. In the next
section we will discuss simulation model for single server, single queue.
monitored. Possible in one time interval, only one customer arrives and only one leaves.
Fixed time increment simulation is generally preferred for continuous simulation. Numerical
methods, where time is taken as independent variable are one such example.
(b) Next Event Increment Simulation: This method is also called Event Oriented
Simulation. In this system, time is incremented when an event occurs. For example in
queuing, when a customer arrives, clock is incremented by his arrival time. In such case
time period for simulation may be stochastic.
Let us consider an example where a factory has a large number of semiautomatic machines. Out
of these machines, none of the machine fails on 50% of the days, whereas on 30% of the days,
one machine fails and on 20% of the days, two can fail. The maintenance staff of the factory, can
repair 65% of these machines in one day depending on the type of fault, 30% in two days and 5%
remaining in three days. We have to simulate the system for 30 days duration and estimate the average
length of queue, average waiting time, and the server loading i.e., the fraction of time for which server
is busy.
The given system is a single server queuing model. Arrival here is failure of the machines and,
while maintenance is the service facility. There is no limit on the number of machines, that is queue
length is infinite. Thus
Expected arrival rate = 0 × 0.5 + 1 × 0.3 + 2 × 0.2 = 0.7 per day
And
Average service time is = 1 × 0.65 + 2 × 0.3 + 3 × 0.05 = 1.4 days
Therefore
Expected service rate = 1/1.4 = 0.714 machines per day
The expected rate of arrival is slightly less than the expected servicing rate and hence the system
can reach a steady state. For the purpose of generating the arrivals per day and the services completed
per day the given discrete distribution will be used.
We generate uniform random numbers between 0 and 1 to generate the failure of machines (arrival
of machines) as under.
0.0 < r ≤ 0.5 arrival = 0 machine
a uniform random number between 0 and 1. Let this number be r = 0.413, which being less than
0.5, from (7.21) we find that no new machine arrives for repair and thus server remains idle.
Day Two: Day is incremented by one, and timer shows 1. By this time server’s idle time is
one day and waiting time is also one day. Again a random number is generated to find new arrival.
Let this number be r = 0.923. From equation (7.21), we find that two machines arrive. Since server
is idle one goes for service and other waits in queue. For simulating service time we again generate
a uniform random number and let this number be r = 0.516. From equation (7.22) we observe that
for this random number repair time (service time) is one day.
Day Three: The situation on third day is: Idle time = 1day, queue = 1, and waiting time = 1
day. We increment day by one and timer becomes 3. Random number of arrival is .571 and thus
one customer comes, and one in queue goes to server. Random number for server is .718 and thus
service time for this machine is 2 days. Today’s customer is in queue. Waiting time becomes 2 days
for this customer who is in queue as server is serving second customer of day 2. Thus waiting time
becomes 2. The process goes on like this.
It is to be noted that when service time is one day, means each machine takes one day for repair
and when service time is 2 days, each machine takes 2 days to be served.
Results of the simulation are given in the Table 7.1.
17 .791 1 2 .145 1 2 27
18 .261 0 1 .801 2 2 28
19 .112 0 1 1 2 29
20 .061 0 0 .081 1 2 29
21 .461 0 0 0 3 29
22 .131 0 0 0 4 29
23 .912 1 1 .161 1 4 30
24 .456 1 1 .881 2 4 31
25 .761 2 2 1 4 33
26 .123 1 1 .531 1 4 34
27 .484 0 0 .981 3 4 34
28 .7101 1 1 4 4 35
29 .533 2 2 1 4 37
30 .901 3 3 .811 2 4 40
Let at time t, (i – 1) customers have arrived into the machine shop and (j – 1) customer have
departed, and i-th customer is about to come and j-th customer is due to depart. If there are some
customers in the queue then
1 ≤ j ≤ i ≤ N
Thus the queue length is (i – j – 1), if i > j. Now next arrival time NAT = CATi and next departure
time (NDT), i.e., the cumulative departure time CDTj of the j-th item is given by,
NDT = CDTj = cumulative arrival time of j + waiting time of j + machining time of j.
= CATj + WTj + STj.
We must now determine which event would take place-whether i would arrive first or j would
depart first. This is decided by comparing NAT with NDT. Now if NAT < NDT, an arrival will
take prior to departure and queue length will increase by one. If NAT > NDT, and (i – j – 1) is
also positive than a departure will take first and queue length will decrease by one. In both these
cases there is no idle time for the server. However if NAT > NDT and the queue length is zero,
then the server will be idle waiting for the i-th item for the duration IDTi = (NAT – NDT). Third
case will be when NAT = NDT, implies then the next arrival and departure take place simulta-
neously and there is no change in the queue length.
This next event simulation procedure for this simple queuing situation is shown in the flow chart
(Fig. 7.3). In the flow chart, the inter arrival times ATi’s and the service times STi’s can be generated
by calling the suitable subroutines. From the inter arrival times, the cumulative arrival times are easily
calculated using the relation,
CATk = CAT k–1 + ATk,
and CAT1 = AT1 = 0
The event times are indicated by the variable time CLOCK.
With this algorithm, one can compute various statistics. This system is demonstrated by following
example.
Example 7.5: The time between mechanic request for tools in a large plant is normally distributed
with mean of 10 minutes and standard deviation of 1.5 minutes. Time to fill requests is also normal
with mean of 9.5 minutes and standard deviation of 1.0 minute. Simulate the system for first 15
requests assuming a single server. Determine the average waiting time of the customers and the
percentage capacity utilization. Develop a computer simulation program and obtain the said results
for 1000 arrivals of customers.
Solution: The normally distributed times, the inter arrival times as well as the service times can
be generated by using normal random number tables (see Appendix-9.1). Normal variable with given
mean (µ) and standard deviation (σ) is given by,
x = µ + σ . (r′)
where r′ is normal random number with mean = 0 and standard deviation unity.
Simulation of the given case is carried out in the Table 7.2. First column gives arrival number.
Second column gives the normal random number generated by Mean Value Theorem, for generating
Inter Arrival Times (IAT). Third column gives inter arrival time, fourth gives cumulative arrival times,
Simulation of Queuing Systems 177
fifth is time for beginning the service, sixth column again is normal random number for service counter,
seventh column is service time, eighth column is service end time, ninth column is customer waiting
time and tenth column is Server Idle time.
For any arrival i, the service begins (SB) either on its arrival cumulative time or the service ending
time of previous arrival, which ever is latest.
SB(i) = Max[SE(i – 1), CAT(i)]
Service ending time is,
SE(i) = SB(i) + ST(i)
Waiting time for customer is,
WT(i) = SE(i – 1) – CAT(i), if possible
Service idle time is,
IT(i) = CAT(i) – SE(i – 1), if possible
178 System Modeling and Simulation
Program for this simulation problem is given below. At the end of 10 simulations,
Average waiting time = 0.2669
Percentage capacity utilization = 80.88
{
sb=cat;
it=sb-se;
cit=cit+it;
}
outfile<<j<<‘\t’<<setprecision(4)<<x1<<‘\t’<<setprecision(4)<<iat<<‘\t’<<setprecision
(4)<<cat<<‘\t’<< setprecision (4)<<sb<<‘\t’<< setprecision (4)<<x2<<‘\t’<<
setprecision (4)<<st<<‘\t’<< setprecision(4)<<se<<‘\t’<<setprecision(4)
<<wt<<‘\t’<<setprecision (4)<<it<<“\n”;
}
awt=cwt/run;
pcu=(cat–cit)*100./cat;
outfile<<“Average waiting time\n”;
outfile<<awt;
outfile<<“Percentage capacity utilization\n”<<pcu;
}
Output of this program is given below. Ten columns are respectively number of arrival, I,
random number r′, inter arrival time IAT, cumulative arrival time CAT, time at which service begins
SB, random number for service time r′, service time ST, time for service ending SE, waiting
time WT, and idle time IT.
Table 7.2: Output of program
1 2 3 4 5 6 7 8 9 10
i r′ IAT CAT SB r′ ST SE WT IT
1 10.72 10.72 10.72 10.72 7.37 7.37 18.09 0 10.72
2 8.493 8.493 19.21 19.21 10.77 10.77 29.98 0 1.123
3 11.68 11.68 30.9 30.9 9.021 9.021 39.92 0 0.913
4 10.74 10.74 41.63 41.63 9.469 9.469 51.1 0 1.714
5 9.095 9.095 50.73 51.1 9.905 9.905 61.01 0.374 1.714
6 9.972 9.972 60.7 61.01 9.988 9.988 70.99 0.306 1.714
7 11.35 11.35 72.05 72.05 9.161 9.161 81.21 0.306 1.057
Contd...
180 System Modeling and Simulation
The problem can be extended to multiple servers with single queue and is given in next section.
Example 7.6: Simulate an M/M/1/∞ queuing system with mean arrival rate as 10 per hour and
the mean service rate as 15 per hour, for a simulation run of 3 hours. Determine the average customer
waiting time, percentage idle time of the server, maximum length of the queue and average length
of queue.
Solution: The notation M/M/1/∞ is Kendal’s notation of queuing system which means, arrival
and departure is exponential, single server and the capacity of the system is infinite.
Mean arrival rate is = 10 per hrs or 1/6 per minute.
Mean departure rate is = 15 per hrs or 1/4 per minute.
In actual simulation the inter arrival times and service times are generated, using exponential
random numbers. Exponential random number R is given by,
R = –(1/τ) ln (1 – r)
where τ is arrival rate. Below we give a program for generating these numbers.
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>//contains rand() function
#include <math.h>
#include <conio.h>
#include<iomanip.h>
void main(void)
{
/* Single server single queue:
arrival and service times are exponentially distributed.
iat is inter arrival time.
count a and count b are counters for server A and B.
qa(qb) are que length of component A(B).
*/
//M/M/1/infinity queue
double r,iat;
double mue=1/6., lemda=1/5., run=180.;
double
Simulation of Queuing Systems 181
clock=0.00,se=0.00,sb=0.00,nat=0.00,cit=0.00,cwt=0.00,st=0.00,it=0.00,wt=0.00;
int q=0,cq=0,k,count=0,qmax=100;
ofstream outfile(“output.txt”,ios::out);
r=rand()/32768.;
cout<<“rand=”<<r<<endl;
iat=(–1./mue)*log(1–r);
nat = nat+iat;
count=count+1;
// cout<<nat<<count<<endl;
// getch();
while(clock<=run) {
if(q>qmax) qmax=q;
outfile<<setprecision(4)<<clock<<‘\t’<<setprecision(4)<<iat<<‘\t’<<
setprecision(4)<<nat<‘\t’<<setprecision(4)<<q<<‘\t’<<setprecision(4)<<sb<<‘
\t’<<setprecision(4)<<st<<‘\t’<<setprecision(4)<<se<<‘\t’<<setprecision(4)<<it<<‘\
t’<<setprecision(4)<<wt<<‘\t’<<setprecision(4)<<cit<<‘\t’<<setprecision(4)<<cwt<<endl;
if(nat>=se){
if(q>0){
wt=q*(se–clock);
cwt=cwt+wt;
q=q–1;
clock=se;
}
else
{ clock=nat;
r=rand()/32768.;
iat=(–1/mue)*log(1–r);
nat = nat=iat;
count=count+1;
}
sb=clock;
it=clock– se;
cit=cit+it;
r=rand()/32768.;
st=(–1/lemda)*log(1–r);
se=sb+st;
}
else
182 System Modeling and Simulation
{
wt=q*iat;
cwt=cwt+wt;
clock=nat;
q=q + 1;
r=rand()/32768.;
st=(–1/lemda)*log(1–r);
nat=nat+iat;
count=count+1;
}
}
outfile<<“Elapsed time=”<<clock<<“Number of arrivals=”<<count<<endl;
outfile<<“Average waiting time/arrival=”<<cwt/count<<endl;
outfile<<“Average server idle time/arrival=”<<cit*100./clock<<endl;
outfile<<“Qmax=”<<qmax<<endl;
cout<<“\n any ’digit”;
cin >> k;
}
Table 7.3: Output of the program
0.00 0.0075 0.007 0 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.007 4.975 4.975 0 0.0075 1.074 1.082 0.0075 0.00 0.0075 0.00
4.975 9.924 9.924 0 4.975 4.397 9.372 3.893 0.00 3.901 0.00
9.924 3.922 3.922 0 9.924 2.156 12.08 0.552 0.00 4.453 0.00
3.922 3.922 7.844 1 9.924 11.31 12.08 0.552 0.00 4.453 0.00
7.844 3.922 11.77 2 9.924 8.653 12.08 0.552 3.922 4.453 3.922
11.77 3.922 15.69 3 9.924 6.864 12.08 0.552 7.844 4.453 11.77
12.08 3.922 15.69 2 12.08 0.9564 13.04 0.00 0.9428 4.453 12.71
13.04 3.922 15.69 1 13.04 9.792 22.83 0.00 1.913 4.453 14.62
15.69 3.922 19.61 2 13.04 6.198 22.83 0.00 3.922 4.453 18.54
19.61 3.922 23.53 3 13.04 3.603 22.83 0.00 7.844 4.453 26.39
22.83 3.922 23.53 2 22.83 1.812 24.64 0.00 9.657 4.453 36.04
23.53 3.922 27.45 3 22.83 0.075 24.64 0.00 7.844 4.453 43.89
24.64 3.922 27.45 2 24.64 0.473 25.12 0.00 3.327 4.453 47.21
25.12 3.922 27.45 1 25.12 2.266 27.39 0.00 0.9585 4.453 48.17
27.39 3.922 27.45 0 27.39 0.797 28.18 0.00 2.266 4.453 50.44
27.45 3.922 31.38 1 27.39 0.907 28.18 0.00 0.00 4.453 50.44
28.18 3.922 31.38 0 28.18 22.32 50.51 0.00 0.7293 4.453 51.17
Contd...
Simulation of Queuing Systems 183
31.38 3.922 35.3 1 28.18 2.95 50.51 0.00 0.00 4.453 51.17
35.3 3.922 39.22 2 28.18 0.634 50.51 0.00 3.922 4.453 55.09
39.22 3.922 43.14 3 28.18 0.023 50.51 0.00 11.77 4.453 74.7
47.06 3.922 50.99 5 28.18 2.373 50.51 0.00 15.69 4.453 90.39
50.51 3.922 50.99 4 50.51 3.793 54.3 0.00 17.22 4.453 107.6
50.99 3.922 54.91 5 50.51 4.233 54.3 0.00 15.69 4.453 123.3
54.3 3.922 54.91 4 54.3 4.603 58.9 0.00 16.58 4.453 139.9
54.91 3.922 58.83 5 54.3 4.672 58.9 0.00 15.69 4.453 155.6
58.83 3.922 62.75 6 54.3 0.909 58.9 0.00 19.61 4.453 175.2
58.9 3.922 62.75 5 58.9 5.439 64.34 0.00 0.4479 4.453 175.6
62.75 3.922 66.67 6 58.9 2.996 64.34 0.00 19.61 4.453 195.2
64.34 3.922 66.67 5 64.34 2.17 66.51 0.00 9.549 4.453 204.8
66.51 3.922 66.67 4 66.51 0.294 66.81 0.00 10.85 4.453 215.6
66.67 3.922 70.59 5 66.51 4.678 66.81 0.00 15.69 4.453 231.3
66.81 3.922 70.59 4 66.81 7.646 74.45 0.00 0.6668 4.453 232
70.59 3.922 74.52 5 66.81 8.112 74.45 0.00 15.69 4.453 247.7
74.45 3.922 74.52 4 74.45 3.668 78.12 0.00 19.29 4.453 267
74.52 3.922 78.44 5 74.45 1.797 78.12 0.00 15.69 4.453 282.6
78.12 3.922 78.44 4 78.12 10.44 88.56 0.00 18.02 4.453 300.7
78.44 3.922 82.36 5 78.12 6.485 88.56 0.00 15.69 4.453 316.4
82.36 3.922 86.28 6 78.12 15.6 88.56 0.00 19.61 4.453 336
86.28 3.922 90.2 7 78.12 13.00 88.56 0.00 23.53 4.453 359.5
88.56 3.922 90.2 6 88.56 3.875 92.43 0.00 15.91 4.453 375.4
90.2 3.922 94.13 7 88.56 0.768 92.43 0.00 23.53 4.453 398.9
92.43 3.922 94.13 6 92.43 3.10 95.53 0.00 15.59 4.453 414.5
94.13 3.922 98.05 7 92.43 1.341 95.53 0.00 23.53 4.453 438.1
95.53 3.922 98.05 6 95.53 9.91 105.4 0.00 9.836 4.453 447.9
98.05 3.922 102.00 7 95.53 1.176 105.4 0.00 23.53 4.453 471.4
102.00 3.922 105.9 8 95.53 7.562 105.4 0.00 27.45 4.453 498.9
105.4 3.922 105.9 7 105.4 9.278 114.7 0.00 27.77 4.453 526.7
105.9 3.922 109.8 8 105.4 28.67 114.7 0.00 27.45 4.453 554.1
109.8 3.922 113.7 9 105.4 40.0 114.7 0.00 31.38 4.453 585.5
113.7 3.922 117.7 10 105.4 4.727 114.7 0.00 35.3 4.453 620.8
114.7 3.922 117.7 9 114.7 2.491 117.2 0.00 9.834 4.453 630.6
117.2 3.922 117.7 8 117.2 1.548 118.8 0.00 22.42 4.453 653
117.7 3.922 121.6 9 117.2 1.764 118.8 0.00 31.38 4.453 684.4
118.8 3.922 121.6 8 118.8 9.167 127.9 0.00 9.904 4.453 694.3
121.6 3.922 125.5 9 118.8 0.1201 127.9 0.00 31.38 4.453 725.
125.5 3.922 129.4 10 118.8 2.357 127.9 0.00 35.3 4.453 761
127.9 3.922 129.4 9 127.9 0.486 128.4 0.00 24.23 4.453 785.2
128.4 3.922 129.4 8 128.4 5.653 134.1 0.00 4.374 4.453 789.6
Contd...
184 System Modeling and Simulation
129.4 3.922 133.3 9 128.4 0.2893 134.1 0.00 31.38 4.453 821
133.3 3.922 137.3 10 128.4 0.0441 134.1 0.00 35.3 4.453 856.3
134.1 3.922 137.3 9 134.1 12.55 146.6 0.00 7.187 4.453 863.5
137.3 3.922 141.2 10 134.1 1.614 146.6 0.00 35.3 4.453 898.8
141.2 3.922 145.1 11 134.1 1.593 146.6 0.00 39.22 4.453 938
145.1 3.922 149.0 12 134.1 4.432 146.6 0.00 43.14 4.453 981.1
146.6 3.922 149.0 11 146.6 5.875 152.5 0.00 18.06 4.453 999.2
149.0 3.922 153.0 12 146.6 9.088 152.5 0.00 43.14 4.453 1042
152.5 3.922 153.0 11 152.5 6.482 159 0 41.49 4.453 1084
153.0 3.922 156.9 12 152.5 3.317 159 0 43.14 4.453 1127
156.9 3.922 160.8 13 152.5 1.149 159 0 47.06 4.453 1174
159.0 3.922 160.8 12 159.0 6.807 165.8 0.00 27.24 4.453 1201
160.8 3.922 164.7 13 159.0 3.16 165.8 0.00 47.06 4.453 1248
164.7 3.922 168.6 14 159.0 3.062 165.8 0.00 50.99 4.453 1299
165.8 3.922 168.6 13 165.8 14.89 180.7 0.00 14.82 4.453 1314
168.6 3.922 172.6 14 165.8 6.821 180.7 0.00 50.99 4.453 1365
172.6 3.922 176.5 15 165.8 0.573 180.7 0.00 54.91 4.453 1420
176.5 3.922 180.4 16 165.8 4.569 180.7 0.00 58.83 4.453 1479
λ
P1 = P0
µ
2
λ 1 1 λ
P2 = P1 +
2µ 2µ
(µP1 − λP0 ) = 2 µ P0
3
λ 1 1λ
P3 = P2 + ( 2µP2 − λ P1 ) = P0
3µ 3µ 3! µ
Thus generalizing above results, the probability that there are n customers in the system,
For n ≥ s is given by,
n
1 λ
Pn = P , for n = 0, 1, 2,..., s ...(7.23)
n! µ 0
1
Input Output
2
(Customers)
Waiting line Service facility (Customers)
n
∞
1 λ
= ∑ n − s P0
n= s s !s µ
186 System Modeling and Simulation
s ∞
1 λ λ
=
s! µ
∑ sµ P0
n= s
s ∞
1 λ λ
=
s! µ
∑ sµ P0
n= s
...(7.24a)
1
P0 =
λ
s
n= s
1 λ µ
n n−s
∞
λ
1 + ∑ + ∑ sµ
...(7.25)
n ! µ s!
n =1 n = s +1
1
P0 =
1 λ
s n
1 λ
s
∑ + ...(7.26)
n = 0 n ! µ s !(1 − λ / µs ) µ
where n = 0, term in the last summation yields the correct value of 1 because of the convention
that n! = 1 when n = 0.
Average length of the queue (LQ) is,
∞
LQ = ∑ ( n − s) P
n=s
n
= ∑ jP
j =0
s+ j
( λ / µ )s
j
∞
λ
= ∑
j= 0
j
s! sµ P0
( λ / µ )s ρ ∞ d
= P0
s!
∑ d ρ (ρ
j=0
j
) ...(7.27)
which is same as
s
λ
µ ρP0
s !(1 − λ / µs ) 2
Average number of customers in the system are equal to average number of customers in the
queue plus average number of customers in the server i.e.,
Simulation of Queuing Systems 187
λ
LS = LQ + ...(7.28)
µ
WS = LS /λ and W Q = LQ /λ ...(7.29)
where LQ is the average number of customers waiting in the queue to be served.
Similarly probability that waiting time is greater than t is,
(λ / µ ) s P0 [1 − e −µ t ( s −1− λ / µ ) ]
P (T > t ) = e −µ t 1 + ...(7.30)
s !(1 − λ /µs ) ( s − 1 − λ /µ )
Example 7.7: In an office there are two typist in a typing pool. Each typist can type an average
of 6 letters/hour. If letters arrive for typing at the rate of 10 letters/hour, calculate,
(a) What fraction of the time are all the two typists busy?
(b) What is the average number of letters waiting to be typed?
(c) What is the waiting time for a letter to be typed?
Solution:
Here λ = 10 letters/hour
µ = 6 letters/hour
s=2
(a) Here we have to find the probability that there are at least two customers in the system,
that is, we want that P(n ≥ 2). To get this
1
P0 =
1
1 + ( λ /µ ) + 2 ( λ /µ ) (1/(1 − λ / 2µ))
2
1
=
1
1 + ( 10/6 ) + 2 ( 10/6 ) (1/(1 − 10/12))
2
= 0.06688
Therefore
( λ /µ )s P0
P (n ≥ 2) =
s !(1 − λ /µ s )
(1.67) 2 × 0.06688
=
2(1 − 10 /12)
2.7889 × 0.06688
= = 0.7971
0.234
this means 79.71% of the time both the typists are busy.
188 System Modeling and Simulation
Note that the probability of one letter in the system (none waiting to be typed and one being
typed) is
2.7889 × 0.06688
P1 = = 0.7971
0.234
(b) Average number of letters waiting to be typed are given by,
( λ /µ )s +1 P0
Lq =
s . s !(1 − λ / µs)2
λ
L = Lq + = 1.15999 + 1.6667 = 2.8267
µ
Therefore,
W = L /λ = 2.8267/10 = 0.283 hr = 13 minutes
Counter 1 Counter 2
k AT k CATk STk,1 CDTk,1 IDTk,1 ST k,2 CDTk,2 IDT k,2 WT k QLk
1 0 0 10 10 0 – – – 0 0
2 9 9 – – – 7 16 9 0 0
3 4 13 12 25 3 – – – 0 0
4 9 22 – – – 20 42 6 0 0
5 4 26 15 41 1 – – – 1 1
6 7 33 – – – 15 48 0 8 1
As the first customer arrives at bank, he directly goes to counter 1. Let service time of first
customer at the service counter is say 10 minutes. Thus within ten minutes first customer leaves
the first counter. Column four gives STk,1 the service time for k-th customer at column 1. Similarly
column seven gives service time at counter 2. Column five gives the cumulative departure time of
k-th customer from counter 1, denoted by CDTk,1. In the sixth column is the idle time IDTk,1 of counter
1, waiting for k-th customer to arrive. Similarly next three columns are defined for service counter
2. The tenth column gives the waiting time (WTk) for k-th customer in the queue. Last column gives
the queue length (QLk) immediately after the arrival of k-th customer.
At time 9 minutes second customer arrives and directly goes to second counter as first counter
is busy. Let service time of second customer at the second counter is 7 minutes, it leaves the counter
at 16 minutes i.e., CDT2,2 = 16. Tenth column gives idle time for second counter i.e., IDT2,2 = 9,
as the second counter was idle, till the second customer arrived.
Third customer arrives at 13 minute and checks if some counter is free, otherwise he waits in
the queue for counter to be free. This is done by comparing the latest values of the cumulative departure
times CDT1,1 and CDT2,2. The smaller of the two indicates when the next facility will be free. Thus
the next departure time is,
MNDT = min(10,16)
Thus Counter one became free at time 10 and is waiting for the customer, who goes to it. Let
its service time be 12 minutes, therefore departure of this customer from counter 1 will be 25 minutes.
This counter remained idle for 3 minutes and queue length is zero.
Next customer arrives at time 22 minutes i.e.,
CAT4 = 22. Since
CAT4 > MNDT(CDT3,1, CDT2,2)
CAT4 – MNDT = 22 – 16 = 6 = IDT4,2
Cumulative departure time at counter number 2 of customer number 4 is given by,
CDT4,2 = CAT4 + ST4,2
= 22 + 20 = 42
Fifth customer arrives at 26 minutes and cumulative departure time of counter 1 is 25 minutes.
This customer goes to counter one and its idle time is one minute.
Similarly when sixth customer arrives at 33 minute no counter is free. Minimum cumulative departure
time is 41 of counter 1. This customer has to wait for eight minutes and then go to counter 1.
190 System Modeling and Simulation
Example 7.8: In an M/D/2/3 system, the mean arrival time is 3 minutes and servers I and II
take exactly 5 and 7 minutes respectively (lambda 1 and lambda 2 in the program) to serve a customer.
Simulate the system for the first one hour of operation, determine the idle time of servers and waiting
time of customers.
Solution: In an M/D/2/3 system, arrivals are distributed exponentially, while the service times
are deterministic. Digit 2 here means there are two servers and 3 means, total number of customers
in the system can not exceed three. The next customer in this case will be returned without service.
The queue discipline will be taken as FIFO.
The mean inter arrival time = 3 minutes
Which means arrival rate λ = 1/3 minutes
For arrival schedule of customers, the random number x, for exponential distribution is generated as,
x = (–1/λ) . ln(1 – r)
where r is a uniform random number between 0 and 1.
Simulation program is given in the Computer program. It is assumed that at time zero, there is
no customer in the queue. To compute the time of first arrival, we compute an exponential random
number.
Program 7.3: Single Queue Two Servers Simulation
// Single queue two servers simulation
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>//contains rand() function
#include <math.h>
#include <conio.h>
#include<iomanip.h>
void main(void)
{ /* M/D/2/3 queuing system.*/
int k,q=0,qmax=3,count=0,counter;
double r, iat,clock=0., nat=0., wt2=0., wt1=0.,it1=0.,it2=0.,cit1=0.,
cit2=0.;
double mean=3.,lambda1=5.,lambda2=4.,se1=0.,se2=0.,run=150;
ofstream outfile(“output.txt”,ios::out);
outfile<<“\n CLOCK IAT NAT SE1 SE2 QUE COUNT
CIT1CIT2 \n”;
// Generate first arrival
while (clock<=run){
//Check the state of arriv7al and update que
r = rand()/32768.0;
iat=(-mean)*log(1–r);
nat=nat+iat;
se1=lambda1;//Service time taken by first server
counter=1;//First customer has come counter=1
Simulation of Queuing Systems 191
outfile.precision(4);
outfile<<clock<<‘\t’<<iat<<‘\t’<<nat<<‘\t’<<se1<<‘\t’<<se2<<‘\t’<<q<<‘\t’
<<count<<‘\t’<<cit1<<cit2<<endl;
//it1 and it2 are idle times for two servers.
while(clock<=run)
{
if(nat<=se1 && nat<=se2){
clock=nat; q=q+1;
r = rand()/32768.0;
iat=(–mean)*log(1–r);
nat=nat+iat; counter=counter+1;
}
else if(se1<=nat && se1<=se2) clock=se1;
else clock=se2;
if (q>qmax){ count=count+1;
q=q–1;
}
if(q>=1 &&se2<=clock)
{
it2=clock–se2;
cit2=cit2+it2;
se2=clock+lambda2;
q=q–1;
}
if(q==0 && se1<=clock)
{
clock=nat;
it1=clock–se1;
cit1=cit1+it1;
se1=nat+lambda1;
se1=nat+lambda1;
r = rand()/32768.0;
iat=(–mean)*log(1–r);
nat=nat+iat;
192 System Modeling and Simulation
counter=counter+1;
}
if(q==0 && se2<=clock)
{
clock=nat;
it2=clock–se2;
cit2=cit2+it2;
se2=nat+lambda2;
r = rand()/32768.0;
iat=(–mean)*log(1–r);
nat=nat+iat;
counter=counter+1;
}
outfile<<clock<<‘\t’<<iat<<‘\t’<<nat<<‘\t’<<se1<<‘\t’<<se2<<‘
\t’<<q<<‘\t’<<count<<‘\t’<<cit1<<‘\t’<<cit2<<endl;
}
outfile.precision(4);
outfile<<“clock=”<<clock<<“cit1=”<<cit1<<“cit2=”<<cit2<<“counter=”<<counter<<endl;
outfile<<“Queuing system M/D/2/3”<<endl;
outfile<<“Mean of the exponential distribution=”<<mean<<endl;
outfile<<“service time of two servers=”<<lambda1<<‘\t’<<lambda 2<<endl;
outfile<<“Simulation run time=”<<clock<<endl;
outfile<<“Number of customers arrived”<<counter<<endl;
outfile<<“Number of customers returned without service”<<count<<endl;
outfile<<“idle time of serverI\n”<<cit1<<endl;
outfile<<“idle time of server II\n”<<cit2<<endl;
outfile<<“Percentage idle time of serverI\n”<<cit1*100/clock<<endl;
outfile<<“Percentage idle time of serverII\n”<<cit2*100/clock<<endl;
}
cout<<“any number”<<endl;
cin>>k;
}
In the Table 7.5, first customer comes at 0.004 minute and goes directly to server1 (SE1) whose
service time is 5 minutes. In second line, second customer comes after 2.487 minutes, so net arrival
time is 2.491 minutes and he goes to server 2 which is idle (cit2) for 0.004 minutes. Third customer
comes at nat = 3.136. Server 1 and server 2, both are busy and hence he waits and queue length
is 1. Fourth customer arrives at cumulative time 8.098 minutes, thus queue length becomes 2. But
second server becomes free at 4.004 minutes and one customer goes to it at 4.004 minutes and clock
time is set at 4.004 minutes and queue length at this time becomes 1. Next line shows clock time
equal to 5 minutes as second customer in queue goes to first server, which becomes free at 5.00
minutes. Similarly further lines can be explained.
Simulation of Queuing Systems 193
Contd...
194 System Modeling and Simulation
EXERCISE
1. Is Poisson’s arrival pattern for queuing is valid for all types of queues? Explain with an
example.
2. How various arrival patterns are generated for the queues? Explain with examples.
(PTU, 2004)
3. Simulate a queue with single queue, two servers. Make your own assumption about the
arrival patterns of customers. (PTU, 2004)
4. Discuss Kendall’s notation for specifying the characteristics of a queue with an example.
5. Jobs arrive at a machine shop at fixed intervals of one hour. Processing time is approxi-
mately normal and has a mean of 50 minutes per job, and a standard deviation of 5 minutes
per job. Simulate the system for 10 jobs. Determine the idle time of operator and job waiting
time, and waiting time of the job. Assume that the first job arrives at time zero. Use the
fixed time incremental model.
6. Repeat the simulation of problem 4, by employing the next event incremental model. Use
the same string of random numbers as in problem 4, and compare the results.
7. An air force station has a schedule of 15 transport flights leaving per day, each with one
pilot. Three reserve pilots are available to replace any pilot who falls sick. The probability
distribution for the daily number of pilots who fall sick is as follows:
Number of pilots sick 0 1 2 3 4 5
Probability 0.15 0.20 0.25 0.15 0.15 0.10
Use the Monte Carlo simulation to estimate the utilization of reserved pilots. What is the
probability of canceling one flight due to non-availability of pilot? Simulate the system for
20 days, 40 days, and 60 days writing a program in C++.
196 System Modeling and Simulation
8. Lallu and Ramu are the two barbers in a barber shop, they own and operate. They provide
two chairs for customers who are waiting to begin a hair cut, so the number of customers
in the shop varies from 0 to 4. For n = 0, 1, 2, 3, 4, the probability Pn that there are
exactly n customers in the shop is P0 = 1/ 16, P1 = 4 /16, P2 = 6 /16, P3 = 4 /16, P4 =1/16.
(a) Calculate LS, queue length.
(b) Determine the expected number of customers being served.
(c) Given that an average of 4 customers per hour arrive and stay to receive a hair cut,
determine WS and WQ.
9. Explain why the utilization factor ρ for the server in a single-server queuing system must
equal 1 = P0 where P0 is the probability of having 0 customer in the system.
10. The jobs to be performed on a particular machine arrive according to a Poisson input
process with a mean rate of two per hour. Suppose that the machine breaks down and
will require one hour to be repaired. What is the probability that the number of new jobs
that will arrive during the time is (a) 0, (b) 2, and (c) 5 or more.
12345678
12345678
12345678
12345678
12345678
12345678
8
12345678
12345678
12345678
SYSTEM DYNAMICS
Control model of autopilot aircraft discussed in chapter one shows how various parts and their activities
are controlled by gyro, and redirects the aircraft in the desired direction. Can we apply control theory
to every day life? It was seen in the example of autopilot that main concern of the control mechanism
is to control the stability and oscillations of the system.
There are many examples in nature where we can also apply control theory. Like engineering
problems, instability and oscillations also occur in nature. Let us take the example of population
explosion. It is observed that if we do not control the population of human being or even any other
specie, it will grow exponentially. This is one of the fields in nature, where control theory is applied
to study the growth of population. In the field of medicines, multiplication of cancerous cells in
human body is another example. Similarly in market there are oscillations in prices of products,
which effect their supply and production. In Physics, decay of radioactive material can also be
modeled with the help of control theory. Although the precision applicable to engineering problems
can not be attained in such problems yet control theory can suggest changes that will improve the
performance of such systems.
In scientific literature studies connected with industrial problems are called Industrial Dynamics
(Geoffrey Gordon), where as study of urban problems is called Urban Dynamics. Similarly control
of environmental problems is called World Dynamics. In all these systems there is no difference in
the techniques to be used to study the system, therefore it is appropriate to call this field as System
Dynamics.
The principal concern of a system Dynamics study is to understand the forces operating on
a system, in order to determine their influence on the stability or growth of the system (Geoffrey
Gordon). Output of such study may suggests some reorganization, or changes in the policy, that
can solve an existing problem, or guide developments away from potentially dangerous directions.
Unlike engineering problems, in this case system dynamics may not produce certain parameters
to improve the performance of system. But this study definitely helps the system analyst to predict
the scenario so that corrective steps can be taken in time.
198 System Modeling and Simulation
Figure 8.1 shows variation of x vs. t for different values of k. Figure shows, higher is the
value of k, more steep is the rise in x. These types of curves are called exponential growth curves
and equation (8.1a) is a mathematical expression for exponential curve. We now apply this law to
the population of monkeys.
Let us take example of population of monkeys in a city. Let their production period is six months.
Assuming that in this city, there are 500 monkeys i.e., 250 couples at time t = 0. If each couple
produces four offspring’s in time = 1 (six months), then proportionality constant k = 2. In first six
months population becomes 1850 (7.4 times) by this relation and in next six months it becomes 13684,
which is equal to 54.7 times of 250.
System Dynamics 199
We can also express function (eq. 8.1) on a semi lag graph. Equation (8.1) can be written as
ln x = ln x0 + 2t ...(8.1b)
which is equation of straight-line in ‘ln x’ (natural log of x) and t. This means that, if we plot this
equation by taking x on log axis and t on simple axis, we will get output as a straight-line, with slope
equal to 2. What are the dimensions of k? From equation (8.1a) it can be seen that dimensions of
k are nothing but 1/time. Sometimes coefficient k is written as 1/T, that is total period under study.
Thus equation (8.1a) can be written as
x = x0et/T
The constant T is said to be time constant, since it provides the measure of growth of x.
Equation 8.2a is plotted in Figure 8.2. From the figure 8.2 we can see that value of x at any
time t, decreases from its initial value and comes down to zero when t is infinity. Decay rate increases
with k.
dx
= k ( X − x) ...(8.3)
dt
with the condition x = 0 at t = 0. Solution of equation (8.3) is as
x = X(1 – e–kt) ...(8.3a)
Figure 8.3 gives the plot of equation (8.3a) for various values of k. This type of curve is sometimes
referred as modified exponential curve. As it can be seen, maximum slope occurs at the origin and
the slope steadily decreases as time increases. As a result of it, the curve approaches the limit more
slowly, and never actually gets the limit. In marketing terms, the sale rate drops as the market
penetration increases. The constant k plays the same role as the growth rate constant as in Exponential
growth model. As k increases, the sale grows more rapidly. As with the growth model, k is sometimes
expressed as equal to 1/T, in which case it can be interpreted as a time constant.
Example 8.1: A builder observes that the rate at which he can sell the houses, depends directly
upon the number of families who do not have a house. As the number of families without house
diminish, the rate at which he sells the houses drops. How many houses in a year can he sell?
Solution: Let H be the potential number of households and y be the number of families with
dy dy
houses. If is the rate at which he can sell the houses, then is proportional to (H – y), i.e.,
dt dt
dy
= k ( H − y ), y = 0 at t = 0
dt
This is nothing but Modified Growth case and solution is
y = H(1 – e–kt)
where H is the potential market.
Example 8.2: Radioactive disintegration
The rate of disintegration of a radioactive element is independent of the temperature, pressure,
or its state of chemical combination. Each element thus disintegrates at a characteristic rate
independent of all external factors. In a radioactive transformation an atom breaks down to give
one or more new atoms.
If to start with t = 0, the number of atoms of A present is a (say). After time t, x atoms will
have decomposed leaving behind (a – x) atoms. If then in a small time interval dt, dx is the number
dx
of atoms which change, the rate of disintegration can be expressed as
dt
dx
= k (a − x ) ...(8.4)
dt
This is called law of mass action and k is called velocity constant or disintegration constant or
transformation constant. Equation (8.4) on integration, with initial condition, x = 0 at t = 0 gives
a−x
ln = – kt
a
or a – x = ae–kt
If we write y = (a – x) and a = y0, we get
y = y0e–kt
202 System Modeling and Simulation
which is same as equation (8.2a). If T is the time when half of the element has decayed i.e., x = a /2 ,
we get T, as
2.303
T= log 2
k
T = 0.693 / k
This means the disintegration of an element to half of its period T depends only on k and is
independent of amount present at time t = 0.
which is the equation for exponential growth curve with proportionality constant equal to kX. Much
later, when the market is almost saturated, the value of x becomes comparable to X, so that it changes
very little with time. The equation for the logistic curve then takes the approximate form
dx
= kX ( X – x )
dt
which is the differential equation for the modified exponential function with a constant kX.
The true differential equation is nonlinear and can be integrated numerically with the boundary
conditions x = 0, when t = 0. Exact analytic solution is tedious and has been given by (Croxton et al.,
1967). Interested students may see the reference. Apart from market trends, many other systems follow
logistic curve, for example population growth can also follow logistic curve (Forrester JW, 1969).
During initial stages, there may be ample resources for the growth of population but ultimately when
resources reduce to scarcity, rate of population comes down.
The model is also applicable to spread of diseases. Initially it spreads rapidly as many acceptors
are available but slowly people uninfected drop and thus growth rate of disease also decreases.
dy
= k1 ( H – y )
dt
dx
= k2 ( y – x) ...(8.6)
dt
Equation (8.6) means the rate of sale of houses at time t, is proportional to (H – y), which is the
number of people who do not have houses. Similarly second equation of (8.6) says, the rate of sale
of air conditioners is proportional to (y – x) i.e., the number of houses which do not have air conditioners
installed so far. Both these equations are modified exponential growth models. In equation (8.6),
we have not taken into account the air conditioners which become unserviceable and require replace-
ment. This factor can be added by modifying second equation of (8.6) as,
dx
= k 2 ( y – x ) – k3 x ...(8.6a)
dt
Equations (8.6) and (8.6a) can easily be computed numerically, when k1, k2 and k3 are known.
204 System Modeling and Simulation
The constant value of k′ shows that the decomposition of KMnO4 is a unimolecular reaction.
dx
= k (a − x)(a − x)
dt
a 2 kt
x=
1 + akt
If we start with different amounts a and b of the reactants, then according to law of mass action,
dx
= k (a − x)(b − x)
dt
with initial conditions, at t = 0, x = 0, one gets,
x
1 −
a
= exp( kt ( a − b))
x
1 −
b
dy 1
= ( H − y)
dt T1
dx 1
= ( y − x) ...(8.7)
dt T2
This is simplistic form of model. Determination of T1 and T2 is not that straight forward. For
example T1 depends on number of factors viz., economic conditions of the people, cost of land, speed
of housing loans available etc. Similarly, T2 depends on again, how many people owning house, can
afford air conditioners, and also it depends on weather conditions and many other factors. In order
206 System Modeling and Simulation
to model such situations, one has to take all the factors into account. In order to model such situation
correctly, feedback of information regarding market trends, is very much essential. Without this
information, air conditioner vendor, may stock air conditioners based on houses available. If he is
not able to install the air conditioners because of other conditions, he will have to stock them for
longer time and bear losses. This becomes another problem called inventory control. Thus all the
parameters have to be taken into account while constructing such models.
Table 8.1
Solution: If we assume that initial value of host and parasite is x0 = 10000 and y0 = 1000, then
following steps will compute the population, by taking t as one day.
xi +1 = xi + (α xi − kxi yi ) dt
yi +1 = yi + ( kxi yi − δyi ) dt
where i = 0, 1, 2, 3,…
Results of computation is given in Table 8.1.
Example 8.4: Babies are born at the rate of one baby per annum for every 20 adults. After a
delay of 6 years, they reach school age. Their education takes 10 years, after which they are adults.
Adults die after an average age of 50 years. Draw a system dynamic diagram of the population and
program the model, assuming the initial number of babies, school going children and adults are
respectively, 300, 3000, and 100,000.
Solution: Let at time t, there are x babies and y adults in the population. Then after time period
δt, increase in number of babies will be,
z . x
δx = δt − δt
20 6
On dividing by δt and as δt → 0, one gets
dx
= z / 20 − x / 6
dt
By similar logic, if y is the population of school going children at time t, then
dy x y
= −
dt 6 10
and if z is the population of adults at time t, then
dz
= kxy – δy
dt
The above equations can be solved numerically by taking initial values for x, y, and z as 300,
3000, and 100,000. Time period t can be taken as one year. Result of this computation is given below.
EXERCISE
1. In the model of house contractor and air conditioners (section 8.5), assume that average
time to sell a house is 4 months, average time to install an air conditioner is 9 months,
and break down of air conditioner occurs, on average 25 month. Take the initial housing
market to be 1000 houses. Assume defected air conditioners to be replaced. Then in a
span of 5 years, how many houses and air conditioners will be sold.
2. In which type of applications, control models are used.
3. Derive an expression for exponential decay models and give one example where it is used.
4. If at time t = 0, radioactive material in a compound is a and proportionality constant is
k, find an expression for half-life of the radioactive part.
5. With reference to market model, is the modified exponential model realistic? Explain in
details, if not, how this model is modified.
6. Give a mathematical model of logistic curve.
7. For a species of animals in an area, the excess of the births over natural deaths causes a
growth rate of a times the current number N. Competition for food causes deaths from
starvation at rate bN 2. Simulate the population growth assuming a = 0.05, b = 0.00001,
and n = 1000 at time t = 0.
8. A certain radio-element has disintegration constant of 16.5 × 106 sec–1. Calculate its half-
life period and average life.
12345678
12345678
12345678
12345678
12345678
9
12345678
12345678
12345678
In chapter seven, we have studied the application of simulation of modeling in various systems where
queuing is involved. Simulation and modeling has application in almost all the branches of science,
especially, where events are stochastic in nature. Inventory control is one such field and will be studied
in this chapter. Whether it is a manufacturing unit or a sale outlet, one of the pressing problem faced
by the management is the control of inventory. Many companies fail each year due to the lack of
adequate control of inventory in their stores. Whether it is raw material used for manufacturing a
product or products waiting for sale, problem arises when, too few or too many items are stored
in the inventory. If the number of items stored are more than what are required, it is a loss of investment
and wastage of storage space which again results in the loss of investments. In the case of excess
inventory, items may depreciate, deteriorate, or become obsolete. But if less number of items are kept
in store, it can result in the loss of sale or reduction in the rate of production, which ultimately results
in the loss of business. In this case there will be loss of profit because of loss of sale and loss of
goodwill due to unfilled demand. Also stock has to be replenished frequently which involves replen-
ishment cost. Then question arises, how much to store in the inventory at a given span of time. This
in turn depends on what is the annual demand and how much time it takes to replenish the inventory
by the supplier. There can be uncertainties, such as strike, weather calamities, price hikes and so
on, in replenishment of inventory. While computing the inventory for storage, all these factors are
to be taken into account. Thus basic problem in inventory control is to optimize the sum of the costs
involved in maintaining an inventory. Backordering is the case when inventory goes to zero and orders
are received for the sale of item, or raw material is required for production. When fresh inventory
arrives, first back orders are completed and then regular orders are entertained. In this case, raw
material is required for production purpose and inventory goes to zero, production stops, which is
a loss to firm. But if inventory is of items, which are for the sale, goodwill of customer is lost in
this case and even there is a possibility that customer will also be lost. Thus backordering case should
always be avoided. In the next section we will discuss basic laws of inventory control models and
also make simulation models for some case studies.
The mathematical inventory models used with this approach can be divided into two broad categories—
deterministic models and stochastic models, according to the predictability of demand involved. If the
210 System Modeling and Simulation
demand in future period can be predicted with considerable precision, it is reasonable to use an
inventory policy that assumes that all forecastes will always be completely accurate. This is the case
of known demand where a deterministic inventory model would be used. However when demand
can not be predicted very well, then in this case stochastic inventory model will be required to be
used, where demand in any period is a random variable rather than a known constant. These cases
are discussed in following sections.
d 2T
and >0 ...(9.4)
dQ 2
First of equation (9.4) gives,
ac0
− + h/2 = 0 ...(9.5)
Q2
or Q = 2ac0 / h
Second of equation (9.4) gives
d 2T
= 2ac0 / Q3 > 0
dQ 2
Thus value of Q given in equation (9.5) gives the minimum value of T.
Equation (9.5) is the well known EOQ formula [21]. Figure 9.1 gives the graphic representation
of inventory depletion and replenishment. The corresponding cycle time t, is given by,
Fig. 9.1: Diagram of inventory level as a function of time for EOQ model.
2c0 a
Q = . ...(9.7)
h (1– a / A)
which is the expression for optimum quantity to be ordered, when arrival rate and depletion rate are
known.
Example 9.3: A small manufacturing company specialises in the production of sleeping bags. Based
on the past records, it is estimated that the company will be able to produce 5000 bags during the
next year if the raw materials is available, when needed. Raw material for each bag costs Rs.50.00.
214 System Modeling and Simulation
Assuming that bags are produced at constant rate during the year of 300 working days, it is estimated
that the annual holding cost of the inventory of raw material is 20% of the raw material cost. Also
each time order is placed, company has to pay Rs.25.00 as reordering cost. If lead time is 7 days,
calculate total annual inventory cost, and total cost.
Solution: Since manufacturer is planning to manufacture 5000 bags in next year, this means
sale (depletion) rate per year is 5000. If an year is taken as time unit, then a = 5000, c0= 25.00,
h = 0.2 × 50 = 10.00, therefore
a Q
TAIC(Q) = c0 . + h.
Q 2
Now Q = 2ac0 / h
Therefore,
a Q a 2ac0 / h
TAIC (Q) = c0 . + h . = c0 +h
Q 2 2ac0 / h 2
2a . c0 . h
= c0 . h. a / 2 +
4
= 25 × 10 × 5000 / 2 + 50 × 10 × 5000 / 4
25 × 5000 10 × 2500
= 5000 × (47.5) + +
2500 2
= 237500.00 + 50.00 + 12500.00 = Rs.250050.00
This cost is less than TC in first case, which is Rs. 251581.14, hence proposal is acceptable.
Example 9.5: In example 9.3, if we put one more condition i.e., arrival rate is 30/day then,
A = 30/day
and a = 5000/300 = 17/day.
2 × 25 × 5000
Then Q=
10 × (1 − 17 / 30)
250000
= = 56818.18 = 238.36
4.4
Suppose we accept the policy of 238 bags each time an order is placed then
hQ a
TC = 5000 × (50) + c0 a /Q + 1 –
2 A
= Rs.250000 + 25 × 50000/238 + 238 × 0.4333
= Rs.250000 + 525.21 + 515.627
= Rs.251040.83
which is slightly less than the total cost in example 9.3, which is 251581.14, hence is economical.
Then let
Q = ordered quantity,
D = annual demand for the product,
B = number of backorders allowed before replenishing inventory,
c0 = ordering cost/order,
h = annual inventory holding cost/unit/time,
p = annual backorder cost/unit/time,
t1 = time for the receipt of an order until the inventory level is again zero,
t2 = time from a zero inventory level until a new order is received,
t3 = time between consecutive orders,
N = D/Q = number of orders/year,
Q – B = inventory level just after a batch of Q units is added.
Since a is the rate of depletion of inventory, during each cycle, the inventory level is positive
for a time (Q – B)/a. The average inventory level during this time is (Q – B)/2 units, and corresponding
cost is h(Q – B)/2 per unit time.
h ( Q – B ) . (Q – B ) h(Q – B) 2
Hence holding cost per cycle = = ...(9.8)
2 a 2a
Similarly shortage occurs for a time B/a. The average amount of shortage during this time is
(0 + B)/2 units, and the corresponding cost is pB/2 per unit time, where p is the cost per unit short
inventory per unit time short.
pB . B p( B)2
Hence shortage cost per cycle = = ...(9.9)
2 a 2a
Reordering cost/cycle = c0, and cost of inventory is cQ. Thus total cost/cycle (TC) is sum of
these costs i.e.,
h(Q – B) 2 pB 2
TC = c0 + + + cQ ...(9.10)
2a 2a
And the total cost per unit time is
c0 a h(Q − B) 2 pB 2
T = + + + ca ...(9.11)
Q 2Q 2Q
Figure 9.2 is the graphic representation of the model.
In this model, there are two decision variables (B and Q), so in order to find optimal values of
∂T ∂T
Q and B, we set partial derivatives and to zero, and solve for Q and B. We get after
∂Q ∂B
differentiation,
∂T − h(Q − B) pB
= + =0
∂B Q Q
h
⇒ B = Q
p+h
∂T − c0 a h(Q − B) h(Q − B) 2 pB 2
and = + − − =0
∂Q Q2 Q 2Q 2 2
2hp
2 p
2
h
2
⇒ 2c0 a = Q − h − p
p + h p + h p + h
2 p p h2
⇒ 2c0 a = Q 2 h − h p + h − p + h
p + h
p+h 2
p h2
⇒ 2 c0 a = Q 2h − h −
p p + h p + h
2ac0 h+ p
or Q=
h p
2ac0 h
and B= ...(9.12)
p p+h
which is optimum values of Q and B.
Example 9.6: Suppose a retailer has the following information available:
a = 350 units/year
c0 = Rs.50 per order
h = Rs.13.75 per unit/time
p = Rs.25 per unit/time
LT = 5 days
To minimize the total annual inventory cost when backordering is allowed, how many units should
be ordered each time an order is placed, and how many backorders should be allowed?
Solution: In this case, optimum order and back order are,
13.75
B= × 63 = 22 units
25 + 13.75
218 System Modeling and Simulation
Thus, the optimal policy is to allow approximately 22 backorders before replenishing the inventory
with approximately 63 units.
t3 + t4 . B
ABC(Q) = c2 ...(9.15)
t5 2
In equation (9.14), (t1 + t2)/t5 is the ratio of time, for which inventory remains positive. Now t5
is the total time in which inventory Q has arrived and depleted. Since arrival rate A is greater than
the depletion rate a, t5 = Q/A. Similarly from Figure 9.3, (t1 + t2) time is the time for which inventory
remained positive, while items were arriving and were being consumed simultaneously. This means,
Q − t1a − B
t1 + t2 =
( A − a)
Inventory Control Models 219
Now t1 is the time in which full EOQ i.e., Q will arrive i.e., t1 = Q/A. Substituting expressions
for t1, t2, t3, t4, t5 in equations (9.13), (9.14) and (9.15) and adding the three costs, one gets
2
c0 D hA . Q A − a − B + pAB 2
TAIC(Q, B) = + ...(9.16)
Q 2Q( A − a ) A 2Q ( A − a )
To find optimum values of Q and B for minimum total annual inventory cost, we differentiate
equation (9.16) with respect to Q and B. After differentiation we get after some algebraic computation,
2co a . h + p
Q=
a ...(9.17)
h(1 − ) p
A
h a
and B= 1 − Q ...(9.18)
h+ p A
variable with a distribution function that is known or can be approximated. The problem is to determine,
how much of the product to have on hand at the beginning of the period to minimize the sum of the
• Cost to purchase or produce enough of the product to bring the inventory upto a certain
level
• Cost of stock-outs (unfilled demands)
• Cost of holding ending inventory
• Since the demand for the product is a random variable, the number of stockouts encountered
and the number of units in ending inventory are also random variables, since they are both
function of demand. Hence, the total inventory cost associated with starting the period with
a given inventory level is a random variable. In this light, we can only hope to determine
the starting inventory level that will minimize the expected value of the three costs that
make up the total inventory cost.
Consider the notation to be used in this section:
X = demand for the product during the given period,
f (x) = distribution function of demand,
F(x) = cumulative distribution function of demand,
Q = order quantity,
c0 = ordering cost per order,
c1 = inventory holding cost/unit of ending inventory,
c3 = cost per item purchased or produced,
c4 = stockout cost/item out of stock,
DIL = desired inventory level at the start of the period,
IOH = initial inventory on hand before placing the order,
TIC(X, DIL) = total inventory cost as a function of demand and the desired inventory level,
TIC(X, DIL, c0) = total inventory cost as a function of demand and the desired inventory level, and set-
up cost c0.
In the next section, we will take the expected value of the total inventory cost with respect to
the demand X, and then determine the value of DIL that will minimize the expected total inventory
cost. We will also assume that backordering is permitted, the delivery rate is finite, and lead time is
zero.
∞ DIL
+ c4
DIL
∫ [( x − DIL)] f ( x) dx ...(9.20)
DIL
+ c4 ∫ [( x − DIL)] f ( x) dx
DIL
which is nothing but sum of cost of Q items, expected holding cost and expected stock-out cost.
Since the expected total inventory cost is not a function of demand, one can take its derivative
with respect to DIL, set it equal to zero, and solve for the optimul DIL i.e.,
DIL ∞
dE[TIC( X , DIL)]
= c3 + c1 ∫ f ( x) dx − c4 ∫ f ( x) dx = 0 ...(9.21)
d (DIL) −∞ DIL
Thus,
c3 + c1 F (DIL) − c4 [1 − F (DIL)] = 0
(c1 + c4 ) F (DIL) = c 4 − c 3
c4 − c3
F (DIL) =
(c1 + c4 )
The inventory level that will minimize the expected total inventory cost is the value of DIL such that
c4 − c3
P( X ≤ DIL) = F (DIL) = ...(9.22)
c1 + c4
The equation (9.22) says F(DIL) represents probability of no stock-outs when the given product
is stocked at the optimum DIL. Likewise
c4 − c3
P ( X > DIL) = 1 − F (DIL) = 1 − ...(9.23)
c1 + c4
represents the probability of atleast one stockout (demand X exceeds DIL).
222 System Modeling and Simulation
Example 9.7: An outdoor equipment shop in Shimla is interested in determining how many pairs
of touring skis should be in stock in the beginning of the skiing season. Assume reordering can not
be done because of the long delay in delivery. Last season was a light year, so the store still has
10 pairs of skis on hand. If
(i)The cost of each pair of skis is Rs.60.
(ii)The retail price is Rs.90.
(iii)The inventory holding cost is Rs.10 per year minus the end of season discount price of
Rs.50.
(iv) The stockout cost is Rs.125 per stockout.
(v) The demand can be approximated with a normal random variable that has a mean of 20
and a variance of 25; X~ N(20, 25).
How many pairs of skis should be stocked at the start of the season to minimize the expected
total inventory cost?
Solution: From the given information
c1 = Rs.10 – Rs.50 = –Rs.40.
c3 = Rs.60
c4 = Rs.125.00
We want to calculate the value of DIL such that
c4 − c3 125 − 60
P( X ≤ DIL) = F (DIL) = = = 0.7647
c1 + c4 −40 + 125
DIL DIL
1 2
or F (DIL) = ∫ f ( x) dx + ∫ e −1/ 2[( x − 20) /5] dx = 0.7647
−∞ −∞ 5 2π
That is, value of DIL, such that the area under the normal curve with mean 20 and variance
25 from –∞ to DIL is equal to 0.7647.
In order to compute the value of DIL, we use normal tables N(0, 1), area under the normal curve
with mean equal to zero and σ equal to 1, the value of z is 0.72 (Appendix 9.1).
DIL z = (DIL − 20)/5
F (DIL) = ∫
−∞
N (20, 25) dx = ∫ N (0,1) dz = 0.7647
−∞
Now-a-days, each addition to inventory and each sale of item is recorded in computer. Whenever
inventory level goes below a specified inventory level, order is placed. For this purpose, several
software packages are available, for implementing such system.
A continuous review inventory system for a particular product, normally will be based on two
critical factors.
DIL = R = reorder point i.e., desired inventory level at the start of the period,
Q = order quantity.
For a manufacturer, managing its finished products inventory, the order will be for a production
run of size Q. For a wholesaler or retailer, the order will be purchase order for Q units of the product.
Thus for these situations, inventory policy would be, whenever the inventory level of the product
drops to R units, place an order for Q units to replenish the inventory.
Such a policy is often called, reorder point, order quantity policy, or (R, Q) policy for short.
Consequently, overall model would be referred as (Q, R) model.
2 Ac0 c2 + h
Q=
h c2
where A now is the average demand per unit time. For the present model, value of Q given above
will only be approximated value.
In order to choose reorder point R, it will be assumed that a stock-out will not occur between
the time an order is placed and the order quantity is received. Thus we denote by L, the management’s
desired probability that stock-out will not occur between the time an order is placed and the order
224 System Modeling and Simulation
quantity is received. This assumption involves working with the estimated probability distribution of
the following random variables.
X = demand during the lead time in filling an order. If the probability distribution of X is
a uniform distribution over the interval from a to b, then set
R = a + L(b – a),
Because then
P( X ≤ R) = L
Since the mean of this distribution is
E(X) = (a + b)/2,
The amount of safety stock (the expected inventory level just before the order quantity is received)
provided by the reorder point R is
Safety stock = R − E ( X )
a+b
= a + L(b − a ) −
2
1
= ( L − )(b − a).
2
When the demand distribution is other than a uniform distribution, the procedure for choosing
R is similar.
General procedure for choosing R.
1. Choose L.
2. Solve for R such that P(X ≥ R) = L
For example, suppose that D has a normal distribution with mean µ and variation σ2. Given the
value of L, the table for the normal distribution given in appendix 9.1, then can be used to determine
the value of R. In particular, one just needs to find the value of (c0)1–L in this table and then plug
into following formula to find R.
R = µ + ( c0 )1+ L . σ
The resulting amount of safety stock is
Safety stock = R – µ = (c0)1–L . σ
To illustrate, if L = 0.75, then (c0)1–L = 0.675, so
R = µ + 0.675 . σ
This provides, safety stock equal to 0.675 times standard deviation.
Example 9.8: A television manufacturing company produces its own speakers, which are used
in the production of television sets. The TV sets are assembled on a continuous production line at
a rate of 8000 per month, with one speaker needed per set. The speakers are produced in batches
because they do not warrant setting up a continuous production line, and relatively large quantities
can be produced in a short time. Therefore, the speakers are placed in the inventory until they are
needed for assembly in to TV sets on the production line. The company is interested in determining
when to produce a batch of speakers and how many speakers to produce in each batch. In order
to solve this problem, several costs must be considered:
Inventory Control Models 225
(i) Each time a batch is produced, a set-up cost of Rs.12000 is incurred. This cost includes
the cost of “tooling up”, and administrative costs.
(ii) The unit production cost of a single speaker is Rs.10.00, independent of batch size.
(iii) The holding cost of speakers is Rs.0.30 per speaker per month.
(iv) Shortage cost per speaker per month is Rs.1.10.
Solution: Originally, there was a fixed demand of speakers i.e., A = 8000/month, to be assembled
into TV sets being produced on a production line at this fixed rate. However, sale of TV sets have
been quite variable, so the inventory level of finished sets has fluctuated widely. To reduce inventory
holding costs for finished TV sets, management has decided to adjust the production rate for the sets
on daily basis to better match the output with the incoming orders.
Demand for speaker in such a case is quite variable. There is a lead time of one month between
ordering a production run to produce speakers for assembly in TV’s. The demand for speakers during
this lead time is a random variable X that has a normal distribution with a mean of 8,000 and a standard
deviation of 2,000. To minimize the risk of disrupting the production line producing the TV sets,
management has decided that the safety stocks for speakers should be large enough to avoid a stockout
during this lead time, 95 percent of the time.
To apply the model, the order quantity for each production run of speakers should be,
2 Ac0 c2 + h
Q = h c2
This is the same order quantity that is found by the EOQ model with planned shortages, where
constant demand rate was assumed. Here management has chosen a service level of L = 0.95, so
that the normal table gives (c0)1–L = 1.645. Therefore reorder point will be,
EXERCISE
Develop a computer simulation model (any language) of the system to determine the optimal
number of news paper copies, which should be procured, so that the expected profit is
maximum. (PTU, 2003)
3. Suppose that the demand for a product is 30 units per month and the items are withdrawn
at a constant rate. The setup cost each time a production run is undertaken to replenish
inventory is Rs.15, production cost is Rs.1.00 per item and the inventory holding cost is
Rs.0.30 per item per month.
(a) Assuming shortages are not allowed, determine how often to make a production run
and what size it should be.
(b) If shortages are allowed but cost Rs.3.00 per item per month, determine how often
to make a production run and what size it should be.
4. Assume the data of example 9.6. What value of DIL will let management be at least
95 percent confident that no demands will go unfilled?
DIL
z = (DIL − 20)/ 5
= ∫
−∞
N (0, 1) dz = 0.95
Inventory Control Models 227
APPENDIX 9.1
0 0000 0040 0080 0120 0160 0199 0239 0279 0319 0359
.1 0398 0438 0478 0517 0557 0596 0636 0675 0714 0754
.2 0793 0832 0871 0910 0948 0987 1026 1064 1103 1141
.3 1179 1217 1255 1293 1331 1368 1406 1443 1480 1517
.4 1554 1591 1628 1664 1700 1736 1772 1808 1844 1879
.5 1915 1950 1985 2019 2054 2088 2123 2157 2190 2224
.6 2258 2291 2324 2357 2389 2422 2454 2486 2518 2549
.7 2580 2612 2642 2673 2704 2734 2764 2794 2823 2852
.8 2881 2910 2939 2967 2996 3023 3051 3078 3106 3133
.9 3159 3186 3212 3238 3264 3289 3315 3340 3365 3389
0 3413 3438 3461 3485 3508 3531 3554 3577 3599 3621
.1 3643 3665 3686 3708 3729 3749 3770 3790 3810 3830
.2 3849 3869 3888 3907 3925 3944 3962 3980 3997 4015
.3 4032 4049 4066 4082 4099 4115 4131 4147 4162 4177
.4 4192 4207 4222 4236 4251 4265 4279 4292 4306 4319
.5 4332 4345 4357 4370 4382 4394 4406 4418 4429 4441
.6 4452 4463 4474 4484 449 4505 4515 4525 4535 4545
.7 4554 4564 4573 8582 4591 4599 4608 4616 4625 4633
.8 4641 4649 4656 4664 4671 4678 4686 4693 4699 4706
.9 4713 4719 4726 4732 4738 4744 4750 4756 4761 4767
0 4772 4778 4783 4788 4793 4798 4803 4808 4812 4817
.1 4821 4826 4830 4834 4838 4842 4846 4850 4854 4857
.2 4861 4864 4868 4871 4875 4878 4881 4884 4887 4890
.3 4893 4896 4898 4901 4904 4906 4909 4911 4913 4916
.4 4918 4920 4922 4925 4927 4929 4931 4932 4934 4836
.5 4938 4940 4941 4943 4945 4946 4948 4949 4951 4952
.6 4953 4955 4956 4957 4959 4960 4961 4962 4963 4964
.7 4965 4966 4967 4968 4969 4970 4971 4972 4973 4974
.8 4974 4975 4976 4977 4977 4978 4979 4979 4980 4981
.9 4981 4982 4982 4983 4984 4984 4985 4985 4986 4986
Contd...
228 System Modeling and Simulation
z 0 1 2 3 4 5 6 7 8 9
0 4987 4987 4987 4988 4988 4989 4989 4989 4990 4990
.1 4990 4991 4991 4991 4992 4992 4992 4992 4993 4993
.2 4993 4993 4994 4994 4994 4994 4995 4995 4995 4995
.3 4995 4995 4995 4996 4996 4996 4996 4996 4996 4997
.4 4997 4997 4997 4997 4997 4997 4997 4997 4997 4998
.5 4998 4998 4998 4998 4998 4998 4998 4998 4998 4998
.6 4998 4998 4999 4999 4999 4999 4999 4999 4999 4999
.7 4999 4999 4999 4999 4999 4999 4999 4999 4999 4999
.8 4999 4999 4999 4999 4999 4999 4999 4999 4999 4999
.9 5000 5000 5000 5000 5000 5000 5000 5000 5000 5000
12345678
12345678
12345678
12345678
12345678
10
12345678
12345678
12345678
COST-EFFECTIVENESS MODELS
Study of life time cost is one of the vital factor involved during the procurement and maintenance of
an equipment and for its efficient use during useful life span. In the present chapter, estimation of Life
Cycle Cost of a military aircraft/missile system will be discussed. Due to its vital importance in a country’s
defence, Life Cycle Cost of a military aircraft/missile system assumes great importance during its
procurement, modification or up-gradation. For taking any decision, whether it is procurement or
deployment, this is the major task before any management. If the weapon has low value, its capabilities
can also be low and on the other hand, also it is not true that a costly weapon is always superior.
A manufacturer will try to highlight qualities of his product and hide its shortcomings. Quite often
the performance of a weapon is too much over rated. Then what is the solution? That is why,
cost-effectiveness studies are required. Cost-effectiveness study of weapon systems is one of the very
important fields of systems analyses. Its need arises, when a new weapon is to be procured or two
weapons are to be compared for their effectiveness as far as their cost and performance is concerned.
In the present chapter a case study, where a surface to surface missile vs. deep penetrating aircraft
is compared for their cost effectiveness (performance effectiveness and cost involved in performing
a typical mission) will be taken up. In order to achieve this a typical mission is assigned to both types
of weapons and cost involved in performing the mission and its merits and demerits are analysed.
It is to be noted that cost-effectiveness study is mainly dependent on the data provided by the
designer or manufacturer. If data is biased, the results of the study can also be incorrect. The
mathematical models used for the study have also very significant role, which are to be chosen very
carefully.
The necessary theoretical background needed for this chapter has been elaborately discussed in
the previous chapters. That is the reason this chapter has been kept as the last chapter. The models
developed in chapters four and five will be utilised in this chapter.
LCC estimation [16, 49–51, 73]. Basic structure of LCC model for aircraft will be discussed in
section 10.3 and that for missile will be discussed in section 10.4. Although method of evaluation
of LCC for any weapon is same, yet two separate models (for aircraft and missile) have been provided
to have deeper understanding of the subject. Aircraft as well as Surface-to-Surface (SS) missile
have many common missions i.e., there are various tasks, which can be performed by both. For
example both can be used for bombing the targets at far-off places in enemy territories. Only difference
is that SS missile can be used once only whereas aircraft can be deployed on a mission again and
again. But on the other hand, missile has a low cost as well as low attrition rate whereas aircraft
has high cost as well as high attrition rate. Also in an aircraft mission, risk to human life (pilot) is
involved, which is not there, in case of attack by missile. There are various other factors, apart from
this, which have to be considered in this study, and will be discussed in this chapter.
FLY-AWAY COST
CIVIL PURCHASE COST
Strike-off Wastage (SOW): This refers to the loss of aircraft during training in peace time.
It varies from 1.5 to 3 aircraft per 10,000 flying hours (where active life of an aircraft is taken as 10,000
hours). The cost of entire aircraft is included to account for the strike-off wastage. For instance
investment costs are the outlays required to introduce the missile system into the operational force
and are related to the size of the force. These also are one-time investments. In this production cost
is not added because it is recurring and depends upon number of missile to be produced. Annual
operating costs are the outlays required to keep the system in operation.
10.3.1 Research and Development (RD) Costing for Each System and Subsystem
The missile subsystems can be divided basically into three categories, viz.
(i) Airframe (A)
(ii) Avionics (L) and
(iii) Engine (E)
Under category A, various systems/sub-systems are, airframe, hydraulic accessories, hydraulic
reservoir, air bottles, missile containers, pneumatics, linkages, propelled feed system, warheads and
other miscellaneous items for integration. Under category L comes, gyro with electronics, accelerometers,
sensor cluster/SD elements, accelerometers electronics, servo controller, pump motor package,
batteries, and cable loom. Under category E comes, LP Engine.
Cost-quantity relationships over the entire production range for the above are to be obtained.
Consider for example sake, one sub-system i.e., airframe in category A. Cost-quality relationship
in respect of the following classes are to be derived and based on these, cost under each class is
obtained.
(i) Initial Engineering (Hrs-cost Eqn.) Cost – (IEi)
(ii) Development Support cost (static test vehicles, – (DSi)
mock ups, test parts and the labour and materials
cost in respect of engineering effort estimated as a
function of initial engineering hours)
(iii) Initial Tooling cost (Hrs-cost Eqn.) – (ITi)
(iv) Manufacturing Labour cost – (ILi)
(v) Manufacturing Material cost – (IMi)
(vi) Sustaining and rate tooling cost
(Maintenance and increased production rates) – (ISTi)
(Hrs-cost)
Similar exercise is to be carried out for each of the 20 (assumed) sub-systems and therefore RD
cost is given as,
20
RD = ∑ Si
i =1
20
= ∑ [ IE i + DSi + ITi + ILi + IM i + ISTi ] ...(10.3)
i =1
(ii) Spares
(iii) Stocks (like personnel supplies, facilities maintenance supplies organisational equipment
supplies)
(iv) Personnel training
(v) Initial travel
(vi) Initial transportation
(vii) Miscellaneous
Cost under each category is considered and is added up in total cost. For example under costing
of “facilities” above three sub-category are being considered.
(a) Ground support systems
(b) Civil works and
(c) Other facilities
(a) Ground support systems of a typical missile alongwith the corresponding cost columns
are reflected in Table 10.1.
16
From Table 10.1, we get the total cost for all ground support systems: GS = ∑ GSi
i =1
Let the costs of (b) Civil works and (c) Other facilities be CW and OF.
234 System Modeling and Simulation
Then the total initial investment cost under facilities category will be II2 = GS + CW + OF.
Let the costs corresponding to the other categories under initial investment be represented by II2
to II7 as per the orders mentioned above.
If n1 number of missiles are being catered, then the initial investment cost can be obtained as II :
7
II = ∑ IIi / ni ...(10.4)
i =1
(AO)T ( A)T ( L )T ( E )T
CTm = C0 K 0 + K1 + K2 + K3 + K4 ...(10.7)
(AO)0 ( A)0 ( L) 0 ( E )0
Then the overall cost of the missile at time T is given by the addition of C Tm , and II i.e.,
4
where Ki’ are the relative weight factors such that, ∑ Ki = 1 subscript “0” denotes base year and
i =0
the aircraft to reach the target, pilot has to update its navigation error with reference to wayside points
during day time. Even with on board night vision IR equipment, target acquisition will become
difficult during night attacks because of limited performance. Further, adverse weather also will make
it difficult to acquire the target during day or night.
60000
55000
50000
45000 lrb 25
40000 lrb 35
35000 lrb 40
lrb 50
Cost
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15
Attrition
In fact even 2% attrition rate is considered to be too high. In actual operation, when an aircraft
encounters a heavy ground air defence, dares not to attack and aborts the mission. Since in this study,
it is assumed that the target has to be defeated, attrition rate is taken as a parameter. In Fig. 10.2,
cost of attack by missiles and a typical strike aircraft is shown for 50% damage, for targets located
at 150km versus different aircraft attrition rates for all compatible types of bombs.
We observe from these figures that up to a certain aircraft attrition rate, the cost of attack by
aircraft is less than that for missile whereas it will be costlier than missile after that value. This value
of attrition rate can be read on x-axis for which missile cost curve and aircraft cost curves intersect.
APPENDIX 10.1
1. Apostol, TM.: Mathematical Analysis, p.74, Addison Wesley, Reading, Mass, 1960.
2. Armament Research and Development Establishment, Training at MBB–Technical Report
(System Group), ARDE, Pune, p. 38, ARDE-141(1982).
3. Balaguruswamy, E: Object Oriented Programming with C++, Tata McGraw-Hill Publishing
Co., New Delhi, India.
4. Balaguruswamy, E: Numerical Methods, Tata McGraw-Hill Publishing Co., New Delhi, India.
5. Ball, Robert, E: The Fundamentals of Aircraft Combat Survivability Analyses and Design,
AIAA Education Series, NY.
6. Banks Jerry, John S. Carson II, Berry L. Nelson and David M. Nicol: Discrete-Event System
Simulation, Prentice-Hall of India Private Ltd, New Delhi, 2002.
7. Barton, DK: The Radar Equations (Radar Vol. 2), Artech House, Dedham, Massachusetts,
1974.
8. Birkhoff, DP, et al.: Explosives with Lined Cavities, J. Appl. Phys. 19, p. 563 (1948).
9. Box CEP and ME Muller: A Note on the Generation of Random Normal Deviates, Ann. Math.
Stat., 29, 1958, 610–11.
10. Breaux, HJ and L S, Mohler: A Computational Procedure for a Class of Coverage Problems
for Multiple Shots, Opns. Res. 19, (1971), pp. 636–644.
11. Brookner, E: Radar Technology, Artech House, Dedham, Massachusetts, 1977.
12. Byren S, Gottfried: Elements of Stochastic Process Simulation, Prentice-Hall Inc., 1984.
13. Chou, PC and AK, Hopkin: Dynamic Response of Materials to Intense Impulsive Loading,
Air Force Material Laboratory, Ohio.
14. Courtnery, P R – Green: Ammunition for the Land Battle. Brassey’s [UK], London, pp. 136–142.
15. Croxton, Frederick E, Dudley J Cowden, and Sidney Klein: Applied General Statistics
(3rd Edition), Englewood Cliffs, NJ,: Prentice Hall, Inc, 1967.
16. DAPCA - A Computer Program for Determining Aircraft Development and Production Costs.
RAND Corp. RM-5221-PR Feb 1967.
17. Design Manuals for Impact Damage Tolerant Aircraft Structures, AGARD-AG-238, 1982,
USA (Classified).
18. Dudewicz, EJ and Mishra, SN: Modern Mathematical Statistics, John Wiley & Sons, New
York (1988), p.144.
19. Encyclopaedia Britannica (Micropaedia), Vol. VIII, 1974, p.70.
242 System Modeling and Simulation
20. Engineering Design Handbook, Army Weapon System Analysis, Part-1, AMRDC, DARCON-
P-706-101 (Nov 1977).
21. Erlenkotter, D: Ford Whitman Harris and the Economic Order Quantity Model, Operations
Research, 38, 937–950, 1990.
22. Forrester Jay, W: Urban Dynamics, Cambridge, Mass., MIT Press, 1969.
23. Gillett, Billy E: Introduction to Operations Research, Tata McGraw-Hill Co. Ltd, New Delhi,
1979.
24. Gnedenko, BV: The Theory of Probability, Mir Publishers, 1969.
25. Gordon, Geoffery: System Simulation, Prentice-Hall of India Pvt Ltd, India, 2004.
26. Griliches Zvi: “Distributed Lags: A Survey,”Econometrica, XXXV, No. 1 (1967), 16–19.
27. Grubbs Frank, E: Expected Target Damage for a Salvo of Rounds with Elliptic Normal
Delivery and Damage Functions, Opns. Res. 16, 1021–1026, (1968).
28. Grubbs Frank, E, Approximate Circular and Noncircular Offset Probabilities of Hitting, Opnd
Res.,12, 51–62 (1962).
29. Held, H: Blast Waves in Air, Propellants, Expl and Pyrotech, 8, pp. 1–7 (1983).
30. Hildebrand, FB: Advanced Calculus for Engineers, pp. 160–62, Prentice-Hall, New York, (1949).
31. Hira, DS: System Simulation, S Chand & Company, New Delhi, 2001.
32. Hovanessian, SA., Radar Detection and Tracking Systems, Pertech House, Dedhan, Massa-
chusetts, 1973, pp. 5–32.
33. Kendall E Atkinson: An Introduction to Numerical Analysis, John Wiley & Sons, NY(1978).
34. Kendall, M and Stuart, A: The Advanced Theory of Statistics, p. 276, Charles Griffin and
Co. Ltd. Lond.
35. Kinney, GF and Graham, KJ: Explosive Shocks in Air, Springer Verlag, New York, 1985.
36. Lee, RC, et al.: Guided Weapons, Brassay’s Defence Publishers, New York.
37. Marsaglia, G and T.A. Bray: A Convenient Method for Generating Normal Variables, SIAM
Review, 6, No. 3, pp. 260–64.
38. Military Hand Book, Survivability, Aircraft, Non-nuclear, Vol 1 to 3, Mil-HDBK, p. 336.
39. Morse, PM and GE, Kimball: Methods of Operations Research, MIT, Massachutes and John
Wiley & Sons (1950).
40. Moss, GM, Leeming DW, and CL, Farrar: Military Ballistics, Brassey’s, London (1983).
41. Narsingh Deo: System Simulation with Digital Computer, Prentice-Hall of India Pvt Ltd, New
Delhi, 2003.
42. Natarajan, R and Wadhwa, LK: Assessment of Damage to Circular and Rectangular Targets
by Bomblike and Gunlike Weapons. First Symp on Warhead Tech., Mar 3–4, pp. 307–312,
TBRL, Chandigarh (1983).
43. Naylor, TJ, Balintfy, JL, Burdick, DS, and K Chu, Computer Simulation Techniques, Wiley,
NY, 1966.
44. Nikolsky, SM: A Course of Mathematical Analysis, Vol. I-II, Mir Publishers, (1977), Moscow.
Bibliography 243
45. Pai, Sh. I.: Compressible Fluid Flow, East West Press, India.
46. Pratap, K Mohapatra, P. Mandal and MC Bora,: Introduction to System Dynamics and
Modelling, University Press India Ltd, (1994).
47. Ramanujam, S and Sashidaran B: An Algorithm for Search of a Minimum TOL Strip in a
Runways Hit by Crater Producing Warheads, ARDE Technical Report 738, (Jun 1983).
48. Rama Rao, K, KPS Murty and MR Patkar: Application of Computer Graphics to Performance
Studies of Missile Warhead, Def. Sci. J., 41 (1), 59–67, (1991).
49. RAND Corporation: A Computer Model for Estimating Development and Procurement Costs
of Aircraft (DAPCA III) RAND Corporation R-1854-PR March 1976.
50. Rand Corp: An Appraisal of Models Used in Life Cycle Cost Estimation for USAF Aircraft
Systems. Rand Corp. R-2287-AF Oct 1978.
51. Reymen DP: Aircraft Design: A Conceptual Approach, AIAA Educational series, (1989).
52. Rhaman, A and SH Zaheer: Operations Res. Quart, Vol. 3, (1952).
53. Rubinstein, Reuven, Y: Simulation and the Monte Carlo Method, John Wiley & Sons, New
York.
54. Sanders, DH, A. Franklin Murphy and Robert J. Eng: Statistics: A Fresh Approach, Intn.
Student edn., McGraw-Hill, Kogakusha Ltd, Tokyo, (1980).
55. Scarborough, James B: Numerical Mathematical Analysis, Oxford and IBH Publishing Co.,
New Delhi.
56. Schroeter, G: Probability of Kill and Expected Destroyed Value if the Underlying Distribution
are Rotationally Symmetric. Ops. Res., 24, pp. 586–591, (1976).
57. Sedov, LI, Similarity and Dimensional Methods in Mechanics, Mir Publishers, Moscow, (1982).
58. Seymour Lipscutz: Theory and Problems of Probability, Metric Edition, Schaum’s Outline
Series, NY, (1987).
59. Shanon, RE: System Simulation – The Art and Science, Prentice Hall Inc., (1975).
60. Singh, VP, T Kunjachan and VS Sasi: Estimation of Target Damage due to Submunition Type
Missile Warhead using Simulation Model/Technique, Def. Sci. J., 47, No. 1, pp. 107–117,
(January 1997).
61. Singh, et al.: Propagation of Spherical Shock Waves in Water, Proc. Ind. Acad. Sci., 3, p.169
(1980).
62. Singh, VP and Singh, SK: Interaction of Explosive Shock with Airborne Cylindrical Targets
of Elliptical Cross Section. Def. Sci. J. 40, 3. pp. 231–242 (July 1990).
63. Singh, VP: Damage Mechanism due to Explosive Loaded Shock Waves, CASSA Monograph
No. CASSA/M-2/88, 1988.
64. Singh, SK and Singh, VP: Mathematical Modelling of Damage to Aircraft Skin Panels Subject
to Blast Loading, Def. Sci. J.,41, 3, 305–316. (July 1991).
65. Singh, VP: On under Water Explosions, Workshop on Terminal Ballistics, 23–24 Dec. (1984),
pp. 88–102, IAT, Pune.
244 System Modeling and Simulation
66. Singh, VP and Yuvraj Sin gh: A Model for Estimation of Attrition for Ground Air Defence
System, Def Sci. J, 41 (4), pp. 417–433, (1991).
67. Singh, VP and Yuvraj Singh: Assessment of Vulnerability of LCA Due to Ground Air Defence
System, CASSA Report No. R-93-1, (1993).
68. Singh, VP, Yuvraj Singh and T Kunjachan, A Vulnerability Study on Recce, Observation and
Attack Helicopters in Modern Battle Field Environment, CASSA Report No. R-96-4, (1996).
69. Singh, VP and Yuvraj Singh: Generalised Model for Aircraft Vulnerability by Different Weapon
Systems, Def. Sci. J, 50, No.1, January 2000, pp.13–23.
70. Singh, VP and T Kunjachen: Damage to Targets of Large Dimensions with Weapons of Small
CEP and Lethal Radius., INT. J. Management & Systems, 15, No. 1, (Jan-April 1999).
71. Skolnik, MI: Introduction to Radar Systems, McGraw-Hill Intn. Co. (1981).
72. Spiegel, MR, Schiller, JJ and Srinivasan, RA: Theory and Problems of Probability and
Statistics, Tata McGraw-Hill Co. Ltd, New Delhi.
73. Srinivasan, NK: A CASSA Model for Life Cycle Cost Estimation for Aircraft, Monograph
on System Analysis for Defence, Sept. 97, CASSA Bangalore, pp.182.
74. Stein, JK: Calculus and Analytic Geometry, McGraw-Hill Intern. Edn. (1977), NY.
75. Tate, A and Ballington, EW: The Physics of Deformations and Flow, McGraw-Hill, p. 573.
76. Timoshenko, SP and S. Woinowsky-Keieger: Theory of Plates and Shells, McGraw Hill,
New York, 1959.
77. Tocher, KD: Art of Simulation, English University Press, London (1963).
78. Treften, FB: History of Operations Research, (in Operations Research in Management, Eds.
JF McLoski and F Treften) John Hopkin Univ Press, Baltimore (1954).
79. Weatherburn, CE: A First Course of Mathematical Statistics, University of Cambridge &
ELBS, India (1961).
80. Weapon Planning Directive No. 1/82, Directorate of Operations (Offensive), IAF. (April 1982).
81. Wilfrid J. Dixon and FJ Massey: Introduction to Statistical Analysis, McGraw-Hill Intn Co.,
Tokyo.
82. William C Guenther and PJ. Terragno: A Review of the Literature on a Coverage Problems,
Ann. Maths. Stat., 35. pp. 232–60, (1964).
83. Willams, JD: An Approximation to the Probability Integral, Ann. Math. Stat., 17, 363–369,
(1946).
84. Yadav, HS and Singh, VP: Converging Shock Waves in Metals, PRAMANA, 18, No. 4,
pp. 331–338 (1982).
85. Yogi, AMN and Sharma, SM: Software Model for Simulation of Target Damage by Stick
Bombing, Silver Jublee Monograph on System Analysis for Defence, CASSA, Bangalore,
pp. 143–156, 25–26 Sept 1997.
86. Zaven A. Karian, and Edward J. Dudewicz: Modern Statistical, Systems and GPSS Simulation,
Computer Science Press, WH Freeman and Co., New York.
87. Zeldovich, YB and Raizer, YP: Physics of Shock Waves and High Temperature Phenomena,
AC Press (NY).
INDEX
A chi-square test 90
Activity 2 Computer models 9
Agner Krarup Erlang 159 Continuity 110
Aim points 55
D
Aiming error 56
Deflection Error Probable 59
Air Defence System 131
Deterministic inventory 210
Aircraft
Direct Attack 71
– combat survivability 65, 132
Disintegration constant 201
– survivability 6
Dispersion 56
Airframe 232
Distributed lag models 19
Anti Aircraft Artillery 67
Dynamic simulation 131
Aural (Sound) Signature 68
Dynamic system 2
Autopilot 2
Avionics (L) 232 E
Electronic Counter Measure 66
B Endogenous 2
Back ordering 213
Energy of Fluids 112
Backorder 215 Engine (E) 232
Binomial Entity 2
– distribution 39 Erlang density function 50
– function 45 Essential Elements Analysis 68
Box-Muller Event 2
– Method 103 Exogenous 2
– Transformation 98 Expected Value 33
Exponential
C – decay models 199
Calling source 162 – distribution 48, 163
Carl Friedrich Gauss 52 – growth model 198
Central limit theorem 97 – probability density 164
246 Index
F Poker’s Method 91
Fixed Time Increment 172 Polygon Method 126, 127
Fluid Flow 110 Prithvi 18
Fuse Functioning 136 Probability Distribution Function 37
Probability of hit 72
G Probability theory 5
Gain 69 Proximity fuse 66
Gamma distribution function 49 Proximity Fuzed Shell 138
Pseudo random numbers 82
I
Pursuit-Evasion Problem 118, 119
Industrial dynamics 6, 197
Intersection 29 Q
IR Signature 68 Quantity of a given substance 204
Queue length 161
K
Queuing theory 6, 162
Kendall’s notation 162
Kolmogrove-Smirnov Test 89 R
Radar cross-section 69
L
Radar Signature 68
Large system 2
Radar transmission 69
Life Cycle Cost 229, 230
Radar Warning Receiver 67
M Random Numbers 80, 81, 89
Marsaglia and Bray Method 98 Random variable 29, 33
Mathematical expectation 51 Random Walk Problem 85
Mathematical models 9 Rejection method 87
Mid Square Random Number Generator 84 Runge-Kutta Method 128
Modified exponential curve 201
Monte Carlo 80
S
Sample points 28
N Sample space 28
Next Event Increment Simulation 173 Shock Waves 115, 116
Normal distribution 52 Signal 69
Signal to noise ratio 69
O
Simulation 1, 18, 80
Operations research 3
Single Server Queue 172
P Single Shot Hit Probability 60, 66, 72, 133, 136
Physical models 9 Sonic barrier 115
Poisson’s distribution 44 Standard deviation 35, 57
Poisson’s arrival pattern 163 State of system 161
Index 247