Evaluation of Characteristic Parameters of Dynamic Models

14 Pages·0.372 MB·English

by  Yankov K.

International Conference on
Information Technologies (InfoTech-2011)
15th – 16th September 2011
St. St. Constantine and Elena resort (Varna), Bulgaria

The forum is organized in the frame of
“Dais of the Science of the Technical University-Sofia, 2011”
and unites the joint events:
25th International Conference on
Systems for Automation of Engineering and Research (SAER-2011)
7th International Workshop on
e-Governance and Data Protection (eG&DP-2011)
1st International Workshop on
Human-Computer Interaction and eLearning Systems (HCIeLS 2011)

PROCEEDINGS

Edited by Prof. Dr. Radi Romansky

Sofia, 2011
4 PROCEEDINGS of the Int’l Conference InfoTech-2011

International Program Committee

Luís BARROSO (Portugal) Gwendal LE GRAND (France)
Dencho BATANOV (Cyprus) Karol MATIAŠKO (Slovakia)
Francesco BERGADANO (Italy) Irina NONINSKA (Bulgaria)
Pino CABALLERO-GIL (Spain) Dimitri PERRIN (Ireland)
Ed F. DEPRETTERE (The Netherlands) Angel POPOV (Bulgaria)
Vassil FOURNADJIEV (Bulgaria) Radi ROMANSKY (Bulgaria)
Georgi GAYDADJIEV (The Netherlands) Giancarlo RUFFO (Italy)
Iliya GEORGIEV (USA) Heather RUSKIN (Ireland)
Adam GRZECH (Poland) Radomir STANKOVIĆ (Serbia)
Luis HERNANDEZ-ENCINAS (Spain) Anastassios TAGARIS (Greece)
Ivan JELINEK (Czech Republic) Ivan TASHEV (USA)
Karl O. JONES (UK) Aristotel TENTOV (Macedonia)
Nikola KASABOV (New Zealand) Dimitar TSANEV (Bulgaria)
Nikola KLEM (Serbia) Michael VRAHATIS (Greece)
Oleg KRAVETS (Russia) Vasilios ZORKADIS (Greece)

HCIeLS Program Committee

Angelos ALEXOPOULOS (Ireland) Lilyana NACHEVA-SKOPALIK (Bulgaria)
Isabel AZEVEDO (Portugal) Boris KRUK (Russia)
Ilona BUCHEM (Germany) Elvira POPESCU (Romania)
Linda CASTAÑEDA (Spain) Ana Elena Guerrero ROLDÁN (Spain)
George IVANOV (UK) Paulo SAMPAIO (Portugal)
Tatyana IVANOVA (Bulgaria)
Technical Board:
Malinka IVANOVA (Bulgaria) Anguelina POPOVA (The Netherlands)
Technical University – Sofia University of Utrecht
[email protected] [email protected]

National Organizing Committee

Chairman: Radi ROMANSKY
Members: Angel POPOV, Dimitar TSANEV, Irina NONINSKA, Todor KOBUROV,
Iva NIKOLOVA, Elena PARVANOVA, Dela STOYANOVA
6 PROCEEDINGS of the Int’l Conference InfoTech-2011
C104 Risk Management and Application of HardFibber Process Bus System 109
A. Petrovski*, V. Fustik*, N. Kiteva Rogleva*, G. Leci** (*Macedonia, **Croatia)
C105 Functional Requirements for Electronic Highway and 115
Risk Analysis for Data Management
V. Fustik*, A. Petrovski*, N. Kiteva Rogleva*, G. Leci** (*Macedonia, **Croatia)
C106 Gene Prediction Using the LZW Data Compression Algorithm 121
Valeria Staneva, Galina Momcheva-Gyrdeva (Bulgaria)
Sub-section “C2”: Information Security and Networking
C201 Security Services of Mobile Telemetry Application Protocol 131
Elitsa Gospodinova, Ivaylo Atanasov, Evelina Pencheva (Bulgaria)
C202 ICT Security Risk Management 139
Biljana Bliznakovska (Macedonia)
C203 Adaptive Anomaly-Based Intrusion Detection System 151
Ljupco Vangelski, Ivan Chorbev, Dragan Mihajlov (Macedonia)
C204 Text Data-Hiding Algorithm 161
Nenad O. Vesić, Dušan J. Simjanović (Serbia)
C205 Implementation of Adaptive Mechanism with Aggregation 169
and Fragment Retransmission for 802.11 Wireless Networks
Valentin Hristov (Bulgaria)
C206 Modeling of VoIP Service 175
Irina Noninska (Bulgaria)
Sub-section “C3”: Automation of System Design and System Investigation
C301 Indirect Method for Retrieving, Storing and Exchange of Big Amount 183
of Raw Data in Distributed Software Systems
Ventseslav Shopov, Vanya Markova (Bulgaria)
C302 Simplifying Trigonometric Expressions using Heuristic 191
Todor Markov (Bulgaria)
C303 Boolean Statisfiability Problem 197
Momchil Peychev (Bulgaria)
C304 Protein Sequence Analysis 205
Khubaib Ahmed Qureshi, Faraz Zaidi, Muhammad Yousuf, Saad Riaz (Pakistan)
C305 Reliability Assessment of Safety-Critical Embedded Software Systems 211
Aleksandar Dimov (Bulgaria)
C306 Downscaling Meteosat Second Generation Thermal Infrared Imagery 217
Stavros Kolios, George Georgoulas, Chrisostomos Stylios (Greece)
C307 Evaluation of Characteristic Parameters of Dynamic Models 225
Kaloyan Yankov (Bulgaria)
C308 Model for Management of Recoursconsumption in Multistage Technology 235
with Periodic Processes in Chemical and Biomechanical Industries
Dragomir Dobrudzhaliev , Boyan Ivanov (Bulgaria)
15-16 September 2011, BULGARIA 225
Proceedings of the International Conference on
Information Technologies (InfoTech-2011)
15-16 September 2011, Bulgaria
EVALUATION OF CHARACTERISTIC PARAMETERS
OF DYNAMIC MODELS
Kaloyan Yankov
Medical Faculty, Trakia University,
Armeiska str., 11, Stara Zagora 6000,
[email protected]
Bulgaria
Abstract. This paper is devoted to the study of dynamic processes in biology and medi-
cine. Characteristic parameters are formulated for a given dynamic process described by
a model consisting of differential equations of I and II order. These parameters are
needed to make a comparative analysis of the system response under different stimulus.
The parameter calculation is done using the inverse function of the model. If it is impo-
ssible to obtain the explicit form of the inverse function, it is approximated with a rever-
sible function that helps to find approximate solutions. In the neighborhood of these
solutions a numerical method is applied to reach a final decision with a preset error.
Key words: system identification, numerical method, modeling, inverse function,
muscle contraction

1. INTRODUCTION

Basic approach in the investigation of systems is by excitation using a test
signal and observing the response over a time interval. When the input influence is a
step or impulse force, the reaction can be identified by ordinary differential equations
(ODE) of I and II order. The methods and the software Korelia for detection and
identification of reaction are described in (Yankov, 2006), (Yankov, 2008), (Yankov,
2009), (Yankov, 2010a). Applying the above described methodology process models
are created after influence of different drugs (Tolekova and Yankov, 2006),
(Tolekova and Yankov, 2008), (Yankov and Tolekova, 2010). From a mathematical
point of view the models fully satisfy the requirements for validity and therefore they
can be used to simulate the reactions of the studied systems.
226 PROCEEDINGS of the International Conference InfoTech-2011
Many researchers focus their interest in identifying geometric parameters using
the graph of the system response. These parameters are usually defined by their
ordinate values, and their respective abscissa values are sought after. For example,
the maximum value of the process and the time to reach its decay time, the total
duration of the process, the times for reaching certain amplitudes identified as part of
the maximum amplitude etc. Though without the completeness and integrity of the
model based on algebraic or differential equations, these parameters are convenient
in many aspects of applied research. Depending on whether one works with
experimental data or with already established model of the process, the approach is
different. In most cases where as input data are given values of the function and look
for the argument the task is to find the inverse function of the model. In my previous
works the calculation of these parameters is not discussed, because the aim was to
obtain a mathematical model of the studied process. If necessary, the coordinates of
points of interest are obtained using the GET-operations of software Korelia
(Yankov, 1997), (Yankov, 1999), (Yankov, 2010b) indicating scattered data or parts
of the model graph with interactive tools.
A mathematical model is a description of a system in terms of equations. It is a
generalization, abstraction of the studied processes or objects. It is desirable that
additional characteristics can be obtained from it, not by specific experimental data,
but because they follow some statistical distribution. In addition, a simulation system
is built from models - there are no data on which to determine additional parameters
of the process.
The aim of this work is to formulate specific points in the development of a
dynamic process and offer algorithms for identifying and calculating them using his
model.

2. EXPONENTIAL MODEL

Models described by ODE I order (Fig.1) are trivial to solve. Their monotonous
ensures the existence of inverse functions, uniqueness of the solution and easiness to
obtain in explicit form.

dy(t)
 ry(t)  kU(t)

dt (1)
y(t ) = C
0 0
Where:

r – the rate constant of the process
U(t) – the input step force
к – the proportionality coefficient

Figure 1. Exponential model
15-16 September 2011, BULGARIA 227
The solution of Eq.(1) is:
y(t) C (1er.t)C er.t
 0 (2)
Where: C = (kU)/r is the infinite asymptote.
∞
The inverse function is easy to obtain from Eq.(2):
1  C  y 
t(y)  ln   (3)
r C C
 
0 
Characteristic process parameters are:
 Steady state level Y and settling time t . After time t the process enters and
s s s
remains within a specified tolerance band ε=δ|С - С |, δ<<1 related to the range
∞ 0
of the process |С - С |:
∞ 0
(C Y ) (C C )
 s  0
For Y the following expression is obtained:
s
Y C (1)C (4)
s 0 
By substituting in Eq.(3) an expression for the settling time is obtained:
1  1 
t  ln  (5)
s r 
 Half-time t = y-1(Y /2). The time when the process reaches half of steady state
1/2 s
level. This is an important parameter used to evaluate the absorption of drugs, the
half life of radioactive materials, etc.
 Furthermore the abscissa values and solution for the area under curve (AUC) will
be sought. This feature in recent years has been the subject of interest in medical
and biological sciences. The AUC is interpreted as the total systemic exposure
and is used as a criterion for evaluation of drug absorption, drug bioavailability
and drug concentrations over time and so on. The AUC in interval [t ,t ] is:
1 2
t2  C C er.t2 er.t1
AUC(t ,t )   C (1er.t)C er.t .dt C (t t ) 0  (6)
1 2  0  2 1 r
t1

3. МODEL WITH LIMITED CAPACITY

A model with limited capacity of resources (Fig.2) is described with the
Verhulst - Pearl equation (7).

dy(t)  y(t)
 r 1 y(t)
 
dt  K  (7)
y(t ) = C
0 0
Where:

r – the rate constant
K – the carrying capacity
Figure 2. Мodel with limited capacity
228 PROCEEDINGS of the International Conference InfoTech-2011
The solution y(t) of Eq.(7) is the logistic model:
K
y(t)
K
1 er.t
(8)
C
0
y(t ) = C
0 0
The inverse function can be explicitly expressed:
1 C C 
t(y) ln 0  0  (9)
 
r y K
 
Characteristic process parameters are:
 The steady state level Y : Y = δC + (1-δ).K , δ <<1
s s 0
 The settling time t = t(Y ).
s s
 The half-time t = t(K/2). If the model is the dose-response relationship (the
1/2
abscissa is the dose d of the studied drug and the ordinate is the probability of
response), it is used to analyze the effect of a drug. d is associated with the
1/2
median effective or the lethal dose (ED50 or LD50) (Yankov, 2010c).
 The AUC in interval [t ,t ] is:
1 2
 
 
t2 K K  C  Ker.t2 
AUC(t1,t2)   K dt  Kt2 t1 r ln C0  Ker.t1  (10)
t11 e(r.t)    0 
C
 
0

4. SECOND ORDER ORDINARY DIFFERENTIAL MODEL

The general form of second order differential equation is:
d2y(t) dy(t)
 2 2y(t)  K 2U(t)

dt2 dt u (11)
y(t ),y’(t ) – initial conditions
0 0
Where: U(t) Input impulse force.
ζ Damping ratio.
ω Natural frequency.
K Gain of the system.
u
Possible reactions and their graphs described by this equation are discussed in
(Yankov, 2009). In a study of biological systems, it appears that the most frequent
are overdamped oscillations, when damping ratio ζ > 1. Therefore, researchers are
concerned mainly with the identification and modeling of such processes. Most
studies of such curves are related to the study of muscle contractions. Models based
on different number of parameters are known. Analytical models based on only two
parameters - contraction time and maximal amplitude are discussed in (Milner-
Brown et al., 1973), (Winter, 1979). These parameters are related to the effectiveness
15-16 September 2011, BULGARIA 229
of the process development. The analytical form proposed in (Raikova and Aladjov,
2004), (Raikova et al., 2007) takes into account four parameters. They can be used
for classification of muscle reactions into four distinct types. The same authors
define six parameter analytical functions in (Raikova et al., 2008) for modeling the
shape of the subtracted contractions with the final aim of analyzing the changes in
the parameters of these contractions. The model identification proposed by the
authors, is a function of the selected geometrical parameters. The chosen system
equations are complex and difficult to identify because of the explicit involvement of
all the six parameters.

Figure 3. Characteristic points in overdamped process

Typical for 2nd order ODEs is that even when the analytical solution is known, it
is impossible to find the inverse function in explicit form. Therefore to obtain the
coordinates for selected specific points is necessary to develop efficient numerical
algorithms.
After analysis of published studies of muscle contractions, the following
characteristic parameters and time intervals in the development of an overdamped
process can be formulated (Fig.3).
 Y . The maximal value of the function.
max
 t . The time of the maximum peak Y .
peak max
 t . The earliest time for which the functions is defined. In most cases t =0.
o o
 t . The measurement of intervals or amplitudes often is required to start from
init
a selected point at the beginning of the graph. In particular, this may be t . In
0
general t is given as a part of the range of the signal, e.g. the time at which
init
when the function increases passes the level y = δ .Y :
init Y max
t = y-1(δ .Y ) , δ << 1
init Y max Y
 t =y-1(Y /2) . Half-rise time. The time for which the process reaches one half
hc max
of its maximal value.
 Y . Steady state level. In practice for the determination of Y the range of the
s s
process is used as a reference:
Y = δ.y +(1+δ).Y , δ<<1
s 0 max
230 PROCEEDINGS of the International Conference InfoTech-2011
 t = y-1(Y ). Settling time. The time when the damped oscillation reaches the
s s
value Y .
s
 t – Half-relaxation time. The time when the process decreases to Y /2.
hr max
The interval evaluations of the process can be obtained from the above
parameters.
 T = t - t . Half-rise (half-contraction) time. The time from the start of
hc hc init
the activity to when the function reaches Y /2.
max
 T = t - t Rise (contraction) time. The time from the start of the
c peak init
process to the it’s maximum.
 T = t - t Sustain time. The period of time in which the system
s hr hc
response remains greater than 0.5Y .
max
 T =t - t Relaxation time. It is the time in which the response
r s peak
decreases from Y to Y .
max s
 T = t - t Total duration of the process.
tot s init
 AUC. For muscle contraction axis is the force of contraction F. From
Newton's second law, the general form of impulse-momentum J relation of a
resultant force F over a time interval t to t .is:
1 2
t2
J(t ,t )  Fdt (12)
1 2
t1
I.e. in the study of muscle contraction AUC (t , t ) is an impulse momentum.
1 2
The inverse function y-1(t) in Eq.(11) cannot be obtained in explicit form. All
characteristic points will be determined as a function of identification parameters
starting from Eq.(11). The principle proposed in (Yankov, 2010a) will be followed.
An approximate solution pa for each parameter p will be determined using
i i
geometric considerations. Then a domain Dom(p) is chosen in the neighborhood of
i
this solution. In this domain numerical procedures are applied for reaching a final
decision with a preset error. Convergences and calculation times depend on the range
of Dom(p). Fourth-order Runge-Kutta method is used as a numerical method.
i

4.1. Characteristic points from the increasing branch of the graph

1. (t , Y ). For overdamped process, when the damping ratio is ζ > 1, a
peak max
private analytical solution (Mastascusa, 1989) is:
K   
y(t)  u e1.t e2.t
(13)

2 2 1
   
1  2 1., 2  2 1.
The first derivative of Eq.(13) is used to find the approximate value ta of this
peak
peak:
15-16 September 2011, BULGARIA 231
dy(t) K   
 u 1.e1.t 2.e2.t  0

dt 2 2 1
ln(1/2)
ta  (14)
peak
12
t is searched in the neighborhood of ta . If assume relative deviation δt ,
peak peak p
the domain is:
Dom(t ): [(1-δt ) tа , (1+δt ) tа ] (15)
peak p peak p peak
A numerical procedure is applied in this range to find the value of the point
(t ,Y ). It is a starting point for calculating the other characteristic points.
peak max
2. t . The solution domain in (+Δt) –neighbourhood of t is:
init 0
Dom(ta ): [t , t +Δt] (16)
init 0 0
Parabolic approximation (Fig.4) is applied for items that are part of the upward
graph. It must be convex and its local extremum coincides with the local extremum
of the graph y(t). Let the equation g (t) of a parabola with a vertical axis be:
p

 
(t h)2  4p g (t)k (17)
p
g (t)[y ,Y ], t[t ,t ],
p 0 max 0 peak
where (h,k) is the peak.
Solving Eq.(17) toward t:
 
t  h 4p k  g(t)
Here h=t ; k=Y .
peak max
The unknown parameter is only p.
Figure 4. Approximations of the process

In this instance, the smaller solution for t is needed. The parabola must pass
through the point (y ,t ): g(t )=y(t ). From Eq.(17) for p the following expression is
0 0 0 0
obtained:
 
1 t t 2
p  0 peak (18)

 
4 Y  y
max 0
For a given value of y , the estimated abscissa tа is:
m m
  Y  y
ta  t  t t max m (19)

peak 0 peak
m Y  y
max 0
t is searched in the neighborhood of ta . If assume relative deviation δt , the
m m m
domain is:
Dom(ta ): [(1-δt ).ta , (1+δt ).ta ] (20)
m m m m m