**Computer Simulation of High‐Frequency Electromagnetic Fields**

Andrey D. Grigoriev

[41] Elsayed K. Optimization of the cyclone separator geometry for minimum pressure drop using Co‐Kriging. Powder Technology. 2015;**269**:409‐424. DOI: 10.1016/j.powtec.2015.

[42] Jones DR, Perttunen CD, Stuckman BE. Lipschitzian optimization without the Lipschitz constant. Journal of Optimization Theory and Applications. 1993;**79**(1):157‐181.

[43] Deb K, Pratap A, Agarwal S, Meyarivan T. A fast and elitist multiobjective genetic algo‐ rithm: NSGA‐II. IEEE Transactions on Evolutionary Computation. 2002;**6**(2):182‐197.

[44] Kestilä A, Tikka T, Peitso P, Rantanen J, Näsilä A, Nordling K, Saari H, Vainio R, Janhunen P, Praks J, Hallikainen M. Aalto‐1 nanosatellite–technical description and mission objec‐ tives. Geoscientific Instrumentation, Methods and Data Systems. 2013;**2**(1):121‐130. DOI:

[45] Khurshid O, Tikka T, Praks J, Hallikainen M. Accommodating the plasma brake experi‐ ment on‐board the Aalto‐1 satellite. Proceedings of the Estonian Academy of Sciences.

[46] Crombecq K, Laermans E, Dhaene T. Efficient space‐filling and non‐collapsing sequen‐ tial design strategies for simulation‐based modeling. European Journal of Operational

09.003

192 Computer Simulation

DOI: 10.1109/4235.996017

10.5194/gi‐2‐121‐2013

2014;**63**(2S):258‐266. DOI: 10.3176/proc.2014.2S.07

Research. 2011;**214**(3):683‐696. DOI: 10.1016/j.ejor.2011.05.032

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/67497

### Abstract

High-frequency and microwave electromagnetic fields are used in billions of various devices and systems. Design of these systems is impossible without detailed analysis of their electromagnetic field. Most of microwave systems are very complex, so analytical solution of the field equations for them is impossible. Therefore, it is necessary to use numerical methods of field simulation. Unfortunately, such complex devices as, for example, modern smartphones cannot be accurately analysed by existing commercial codes. The chapter contains a short review of modern numerical methods for Maxwell's equations solution. Among them, a vector finite element method is the most suitable for simulation of complex devices with hundreds of details of various forms and materials, but electrically not too large. The method is implemented in the computer code radio frequency simulator (RFS). The code has friendly user interface, an advanced mesh generator, efficient solver and post-processor. It solves eigenmode problems, driven waveguide problems, antenna problems, electromagnetic-compatibility problems and others in frequency domain.

Keywords: electromagnetics, numerical methods, computer simulation, microwaves, cellular phones

## 1. Introduction

High-frequency electromagnetic fields are used now in telecommunications and radar systems, astrophysics, plasma heating and diagnostics, biology, medicine, technology and many other applications. Special electromagnetic systems excite and guide these fields with given time and space distribution. A designer or a researcher of such systems ought to know in detail their electromagnetic field characteristics. This goal can be achieved or by experimental study,

© 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

often too long in time and expensive, or by computer simulation. The last choice becomes more and more preferable with fast progress of computational electrodynamics and computer efficiency.

All macroscopic electromagnetic phenomena are governed by Maxwell's equations. Unfortunately, these remarkable equations have so many solutions, that choice of the one satisfying given initial and boundary conditions (BC) often becomes very difficult problem. A number of commercial computer codes based on numerical solution of Maxwell's equations are available at this time. These codes make possible high-frequency electromagnetic field simulation.

Researches have not yet created a universal code, efficiently simulating electromagnetic field excited by an arbitrary technical or nature source in an arbitrary medium. Some codes are more suitable for solving one kind of problems and other codes—another kind. Hence, developing new, more universal and efficient computer codes is an actual task. On the other side, a designer of electromagnetic devices and systems has to choose most efficient computer code for solving his/her particular problem. He/she can do right choice only if he/she understands the basics of a numerical method used in the given code.

The goal of the presented chapter is to formulate electromagnetic problems and to describe in short prevailing numerical methods of its solving. As an example, the chapter also gives more detailed description of the radio frequency simulator (RFS) computer code, developed in collaboration of Saint-Petersburg State Electrotechnical University and LG Russian R&D Centre. Some results obtained by means of this code demonstrate its accuracy and efficiency.

The author hopes that this chapter would be useful for researchers and designers of modern telecommunication devices and systems.

## 2. Basic equations

Maxwell's equations are the basic ones, describing macroscopic electromagnetic fields in an arbitrary medium. In modern notation, these equations have the form

$$
\nabla \times \mathbf{H} = \frac{\partial \mathbf{D}}{\partial t} + \mathbf{J} \tag{1}
$$

D ¼ εE ¼ ε0εrE (5)

Computer Simulation of High‐Frequency Electromagnetic Fields

B ¼ μH ¼ μ0μrH (6)

n · ðE<sup>2</sup> � E1Þ ¼ 0; n · ðH<sup>2</sup> � H1Þ ¼ J<sup>s</sup> (7)

n · E ¼ 0 (8)

n · H ¼ 0 (9)

e<sup>n</sup> · ð∇ · ĖÞ ¼ 0 (10)

ωμ=ð2σ<sup>Þ</sup> <sup>p</sup> .

� �E<sup>τ</sup> <sup>¼</sup> <sup>0</sup> (11)

<sup>=</sup>ð4πc<sup>2</sup>Þ, <sup>μ</sup><sup>0</sup> <sup>¼</sup> <sup>4</sup><sup>π</sup> � <sup>10</sup>�<sup>7</sup> are

http://dx.doi.org/10.5772/67497

195

Here, <sup>ε</sup>, <sup>μ</sup> are the absolute permittivity and permeability, <sup>ε</sup><sup>0</sup> <sup>¼</sup> 107

The most frequently used boundary conditions (BC) are as follows:

• Similarly, on the surface of perfect magnetic conductor (PMC)

• Approximate Leontovich boundary condition holds on the impedance surface:

k0η0μ<sup>r</sup>

• We define radiation (adsorption) condition on the surface throw where radiation propagates without reflections. Higdon [1] proposed general absorption boundary conditions (ABCs) theory of arbitrary order of approximation. For a wave, propagating under arbi-

> ∂ ∂t þ ξ

where ξ is constant, providing solution stability and u is the phase velocity of the wave. Higher

<sup>e</sup><sup>n</sup> · <sup>ð</sup>e<sup>n</sup> · <sup>Ė</sup>Þ � <sup>j</sup> Zs

where Zs is the surface impedance. For metals, Zs ¼ ð<sup>1</sup> <sup>þ</sup> <sup>j</sup><sup>Þ</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

trary angle ϕ to the border x ¼ Const, the first-order Higton's ABC are

order ABC can be constructed by multiplying several first-order operators (11).

cosϕ u

∂ ∂x þ

the initial moment of time.

• On the surface separating two dielectrics:

J<sup>s</sup> is the surface electric current density.

• On the surface of perfect electric conductor (PEC):

the dielectric and magnetic constants (we use the SI units in this chapter) and εr, μ<sup>r</sup> are the relative permittivity and permeability, which can be scalars or tensors, depending on medium properties. Equations (1), (2), (5) and (6) form a system of 12 scalar differential equations of the first order with 12 unknowns—components of E, H, D,B vectors, which are functions of space coordinates and time. To solve this system, one needs to define initial and boundary conditions. These conditions together with the equations form an electrodynamics problem (EMP). Initial conditions define electric and magnetic field intensities in the computational region at

where n is the unit normal to the surface directed from the first medium to the second and

$$
\nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} \tag{2}
$$

$$\nabla \cdot \mathbf{D} = \rho \tag{3}$$

$$\nabla \cdot \mathbf{B} = 0 \tag{4}$$

In these formulas, J, ρ are the electric current and electric charge densities (field sources), E, H are the electric and magnetic field intensities (or simply electric and magnetic fields), D,B are the electric and magnetic flux densities, ∇ is Hamilton's differential operator and · is the sign of vector and scalar products.

The constitutive relations couple flux densities and field intensities:

$$\mathbf{D} = \varepsilon \mathbf{E} = \varepsilon\_0 \varepsilon\_r \mathbf{E} \tag{5}$$

$$\mathbf{B} = \mu \mathbf{H} = \mu\_0 \mu\_r \mathbf{H} \tag{6}$$

Here, <sup>ε</sup>, <sup>μ</sup> are the absolute permittivity and permeability, <sup>ε</sup><sup>0</sup> <sup>¼</sup> 107 <sup>=</sup>ð4πc<sup>2</sup>Þ, <sup>μ</sup><sup>0</sup> <sup>¼</sup> <sup>4</sup><sup>π</sup> � <sup>10</sup>�<sup>7</sup> are the dielectric and magnetic constants (we use the SI units in this chapter) and εr, μ<sup>r</sup> are the relative permittivity and permeability, which can be scalars or tensors, depending on medium properties. Equations (1), (2), (5) and (6) form a system of 12 scalar differential equations of the first order with 12 unknowns—components of E, H, D,B vectors, which are functions of space coordinates and time. To solve this system, one needs to define initial and boundary conditions. These conditions together with the equations form an electrodynamics problem (EMP).

Initial conditions define electric and magnetic field intensities in the computational region at the initial moment of time.

The most frequently used boundary conditions (BC) are as follows:

• On the surface separating two dielectrics:

often too long in time and expensive, or by computer simulation. The last choice becomes more and more preferable with fast progress of computational electrodynamics and computer effi-

All macroscopic electromagnetic phenomena are governed by Maxwell's equations. Unfortunately, these remarkable equations have so many solutions, that choice of the one satisfying given initial and boundary conditions (BC) often becomes very difficult problem. A number of commercial computer codes based on numerical solution of Maxwell's equations are available at this time. These codes make possible high-frequency electromagnetic field simulation.

Researches have not yet created a universal code, efficiently simulating electromagnetic field excited by an arbitrary technical or nature source in an arbitrary medium. Some codes are more suitable for solving one kind of problems and other codes—another kind. Hence, developing new, more universal and efficient computer codes is an actual task. On the other side, a designer of electromagnetic devices and systems has to choose most efficient computer code for solving his/her particular problem. He/she can do right choice only if he/she understands

The goal of the presented chapter is to formulate electromagnetic problems and to describe in short prevailing numerical methods of its solving. As an example, the chapter also gives more detailed description of the radio frequency simulator (RFS) computer code, developed in collaboration of Saint-Petersburg State Electrotechnical University and LG Russian R&D Centre. Some results obtained by means of this code demonstrate its accuracy and efficiency.

The author hopes that this chapter would be useful for researchers and designers of modern

Maxwell's equations are the basic ones, describing macroscopic electromagnetic fields in an

<sup>∇</sup> · <sup>H</sup> <sup>¼</sup> <sup>∂</sup><sup>D</sup>

<sup>∇</sup> · <sup>E</sup> ¼ � <sup>∂</sup><sup>B</sup>

In these formulas, J, ρ are the electric current and electric charge densities (field sources), E, H are the electric and magnetic field intensities (or simply electric and magnetic fields), D,B are the electric and magnetic flux densities, ∇ is Hamilton's differential operator and · is the sign

<sup>∂</sup><sup>t</sup> <sup>þ</sup> <sup>J</sup> (1)

<sup>∂</sup><sup>t</sup> (2)

∇ � D ¼ ρ (3)

∇ � B ¼ 0 (4)

arbitrary medium. In modern notation, these equations have the form

The constitutive relations couple flux densities and field intensities:

the basics of a numerical method used in the given code.

telecommunication devices and systems.

2. Basic equations

of vector and scalar products.

ciency.

194 Computer Simulation

$$\mathbf{n} \times (\mathbf{E}\_2 - \mathbf{E}\_1) = 0; \quad \mathbf{n} \times (\mathbf{H}\_2 - \mathbf{H}\_1) = \mathbf{J}\_s \tag{7}$$

where n is the unit normal to the surface directed from the first medium to the second and J<sup>s</sup> is the surface electric current density.

• On the surface of perfect electric conductor (PEC):

$$\mathbf{n} \times \mathbf{E} = 0 \tag{8}$$

• Similarly, on the surface of perfect magnetic conductor (PMC)

$$\mathbf{n} \times \mathbf{H} = 0\tag{9}$$

• Approximate Leontovich boundary condition holds on the impedance surface:

$$\mathbf{e}\_{\boldsymbol{\pi}} \times (\mathbf{e}\_{\boldsymbol{\pi}} \times \dot{\mathbf{E}}) - j \frac{Z\_s}{k\_0 \eta\_0 \mu\_r} \mathbf{e}\_{\boldsymbol{\pi}} \times (\nabla \times \dot{\mathbf{E}}) = \mathbf{0} \tag{10}$$

where Zs is the surface impedance. For metals, Zs ¼ ð<sup>1</sup> <sup>þ</sup> <sup>j</sup><sup>Þ</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ωμ=ð2σ<sup>Þ</sup> <sup>p</sup> .

• We define radiation (adsorption) condition on the surface throw where radiation propagates without reflections. Higdon [1] proposed general absorption boundary conditions (ABCs) theory of arbitrary order of approximation. For a wave, propagating under arbitrary angle ϕ to the border x ¼ Const, the first-order Higton's ABC are

$$\mathbf{E}\left(\frac{\partial}{\partial x} + \frac{\cos\varphi}{u}\frac{\partial}{\partial t} + \xi\right)\mathbf{E}\_t = \mathbf{0} \tag{11}$$

where ξ is constant, providing solution stability and u is the phase velocity of the wave. Higher order ABC can be constructed by multiplying several first-order operators (11).

Another expression for ABC can be derived from the wave equation. For the wave propagating along the normal e<sup>n</sup> to the border, we get

$$(\nabla \times \dot{\mathbf{E}})\_{\tau} + jk\_0 \sqrt{\varepsilon\_r \mu\_r} \text{ e}\_{n \times} \times \dot{\mathbf{E}} = \mathbf{0} \tag{12}$$

2. Analytical treatment—formulating of equations, initial and boundary conditions, geometrical description of the computational region and filling medium properties. The choice of the numerical solution method, transforming equations to the form, most suitable for the

Computer Simulation of High‐Frequency Electromagnetic Fields

http://dx.doi.org/10.5772/67497

197

3. Problem discretization—transfer from continuous functions to discrete ones and from functional equations to the system of linear algebraic equations (SLAE), in the certain

4. Algebraic solution—choosing the most efficient numerical method and solving the SLAE

5. Post-processing—calculation of fields, characteristics and parameters of the electromag-

Each stage of the solution adds its own contribution to the total solution error. The first step adds the so-called inherent error arising due to inaccuracy of input data. This error cannot be eliminated on the next stages of solution. The second stage adds mathematical model error caused by imperfect adequateness of the model to the real physical process. Problem discretization adds the so-called numerical method error, the value of which depends on the quality of the discretization process. At last, computational error arises on stages 4 and 5 due to finite accuracy of numbers presentation in a computer and finite number of operations. With progress in the computational mathematics and computer, the main error sources move from the uncontrolled first stages to the

The solution of real-life electromagnetic problem is a very complicated task. There are no universal methods, capable to solve efficiently an arbitrary problem. Hence, a number of numerical methods were elaborated, each of them is the most efficient for its particular range

The first group solves problems involving Maxwell's equation (1)–(4), containing time as an independent variable. This group makes possible to find the solution in time domain (TD) with arbitrary time dependence of fields. This group of methods is the most suitable for solving non-linear problems. If the problem is linear, Fourier transform can be used to find frequency

The most popular method of this group is the Finite Difference Time Domain (FDTD) method, proposed by Yee [4, 5]. Another method, successfully used in time domain, is the Finite Integration Technique (FIT), firstly proposed by Weiland [6]. The Transmission Line Matrix (TLM) method, proposed by Johns and Beurle [7], also works efficiently in time domain. Its detailed description can be found in [8]. These methods were implemented in a number of computer codes, such as SPEAG SEMCAD™ (FDTD), CST Studio Suite™ (FIT, TLM) and

fourth and fifth stages, where we can often predict the error value a-priory.

of problems. The numerical methods divide on two large groups.

chosen method, a-priory analysis of equations and their solutions properties.

sense approaching the initial problem.

netic system and visualization of the results.

4. Numerical methods classification

spectrum of the solution.

many others.

with the prescribed accuracy.

There are two types of EMPs: an inner problem, when the solution is defined in a closed space region with certain boundary conditions on its border, and an outer problem presuming existence of the solution in the unbounded space excluding some regions with prescribed conditions on their boundaries. We ought also to define field sources in the computational region and initial conditions—electrical and magnetic fields at some moment of time. The proper initial and boundary conditions guarantee existence and uniqueness of the solution [2, 3].

Equations (1)–(4) presume arbitrary time dependence of sources and fields. But very often field dependence on time expresses by the harmonic law: a ¼ am cos ðωt þ ϕÞ, where a is any component of the field, ω ¼ 2πf is the angular frequency, f is frequency and ϕ is the initial angle. In this case, we can simplify Maxwell's equations, using notation

$$\mathbf{E}(\mathbf{r},t) = \text{Re}[\dot{\mathbf{E}}(\mathbf{r})\mathbf{e}^{j\omega t}],\tag{13}$$

where <sup>Ė</sup>ðrÞ ¼ Exe<sup>j</sup>ϕ<sup>x</sup> <sup>e</sup><sup>x</sup> <sup>þ</sup> Eye <sup>j</sup>ϕ<sup>y</sup> <sup>e</sup><sup>y</sup> <sup>þ</sup> Eze<sup>j</sup>ϕze<sup>z</sup> is the complex electric field amplitude (phasor). Magnetic field H is presented similarly. Using these notations, we can write Maxwell's equations for phasors:

$$\nabla \times \dot{\mathbf{H}} = j\omega \dot{\epsilon} \dot{\mathbf{E}} + \dot{\mathbf{J}} \tag{14}$$

$$\nabla \times \dot{\mathbf{E}} = -j\omega\dot{\mu}\dot{\mathbf{H}}\tag{15}$$

$$\nabla \cdot (\dot{\varepsilon} \dot{\mathbf{E}}) = \dot{\rho} \tag{16}$$

$$\nabla \cdot (\dot{\mu}\dot{\mathbf{H}}) = \mathbf{0} \tag{17}$$

Here, ε\_ ¼ ε0ε\_r, μ\_ ¼ μ0μ\_ <sup>r</sup> are the complex absolute permittivity and permeability and ε\_<sup>r</sup> ¼ ðε<sup>0</sup> <sup>r</sup> <sup>þ</sup> <sup>j</sup>ε″ <sup>r</sup>Þ, μ\_ <sup>r</sup> ¼ ðμ<sup>0</sup> <sup>r</sup> <sup>þ</sup> <sup>j</sup>μ″ <sup>r</sup>Þ are the relative complex permittivity and permeability. The equation systems (14)–(17) are simpler than the original one because it does not contain time as independent variable and has only six unknown functions.

#### 3. Basis stages of electromagnetic problem solution

Solution of a given electromagnetic problem can be divided on several subsequent steps (stages):

1. Problem formulation—defining the goal of the computation, necessary input and output data, admissible inaccuracy of the results.


Each stage of the solution adds its own contribution to the total solution error. The first step adds the so-called inherent error arising due to inaccuracy of input data. This error cannot be eliminated on the next stages of solution. The second stage adds mathematical model error caused by imperfect adequateness of the model to the real physical process. Problem discretization adds the so-called numerical method error, the value of which depends on the quality of the discretization process. At last, computational error arises on stages 4 and 5 due to finite accuracy of numbers presentation in a computer and finite number of operations. With progress in the computational mathematics and computer, the main error sources move from the uncontrolled first stages to the fourth and fifth stages, where we can often predict the error value a-priory.

## 4. Numerical methods classification

Another expression for ABC can be derived from the wave equation. For the wave propagating

ffiffiffiffiffiffiffiffiffi εrμ<sup>r</sup>

There are two types of EMPs: an inner problem, when the solution is defined in a closed space region with certain boundary conditions on its border, and an outer problem presuming existence of the solution in the unbounded space excluding some regions with prescribed conditions on their boundaries. We ought also to define field sources in the computational region and initial conditions—electrical and magnetic fields at some moment of time. The proper initial and boundary conditions guarantee existence and uniqueness of the solution

Equations (1)–(4) presume arbitrary time dependence of sources and fields. But very often field dependence on time expresses by the harmonic law: a ¼ am cos ðωt þ ϕÞ, where a is any component of the field, ω ¼ 2πf is the angular frequency, f is frequency and ϕ is the initial angle. In

<sup>E</sup>ðr,tÞ ¼ Re½ĖðrÞe<sup>j</sup>ω<sup>t</sup>

Magnetic field H is presented similarly. Using these notations, we can write Maxwell's equa-

∇ · Ė ¼ �jωμ\_ H

Here, ε\_ ¼ ε0ε\_r, μ\_ ¼ μ0μ\_ <sup>r</sup> are the complex absolute permittivity and permeability and

equation systems (14)–(17) are simpler than the original one because it does not contain time as

Solution of a given electromagnetic problem can be divided on several subsequent steps

1. Problem formulation—defining the goal of the computation, necessary input and output

∇ � ðμ\_ H �

<sup>¼</sup> <sup>j</sup>ωε\_ <sup>Ė</sup> <sup>þ</sup> \_

�

∇ · H �

<sup>p</sup> <sup>e</sup><sup>n</sup> · · <sup>Ė</sup> <sup>¼</sup> <sup>0</sup> (12)

�, (13)

J (14)

∇ � ðε\_ ĖÞ ¼ ρ\_ (16)

<sup>r</sup>Þ are the relative complex permittivity and permeability. The

Þ ¼ 0 (17)

(15)

<sup>j</sup>ϕ<sup>y</sup> <sup>e</sup><sup>y</sup> <sup>þ</sup> Eze<sup>j</sup>ϕze<sup>z</sup> is the complex electric field amplitude (phasor).

ð∇ · ĖÞ<sup>τ</sup> þ jk<sup>0</sup>

this case, we can simplify Maxwell's equations, using notation

along the normal e<sup>n</sup> to the border, we get

[2, 3].

196 Computer Simulation

where <sup>Ė</sup>ðrÞ ¼ Exe<sup>j</sup>ϕ<sup>x</sup> <sup>e</sup><sup>x</sup> <sup>þ</sup> Eye

tions for phasors:

ε\_<sup>r</sup> ¼ ðε<sup>0</sup>

(stages):

<sup>r</sup> <sup>þ</sup> <sup>j</sup>ε″

<sup>r</sup>Þ, μ\_ <sup>r</sup> ¼ ðμ<sup>0</sup>

<sup>r</sup> <sup>þ</sup> <sup>j</sup>μ″

data, admissible inaccuracy of the results.

independent variable and has only six unknown functions.

3. Basis stages of electromagnetic problem solution

The solution of real-life electromagnetic problem is a very complicated task. There are no universal methods, capable to solve efficiently an arbitrary problem. Hence, a number of numerical methods were elaborated, each of them is the most efficient for its particular range of problems. The numerical methods divide on two large groups.

The first group solves problems involving Maxwell's equation (1)–(4), containing time as an independent variable. This group makes possible to find the solution in time domain (TD) with arbitrary time dependence of fields. This group of methods is the most suitable for solving non-linear problems. If the problem is linear, Fourier transform can be used to find frequency spectrum of the solution.

The most popular method of this group is the Finite Difference Time Domain (FDTD) method, proposed by Yee [4, 5]. Another method, successfully used in time domain, is the Finite Integration Technique (FIT), firstly proposed by Weiland [6]. The Transmission Line Matrix (TLM) method, proposed by Johns and Beurle [7], also works efficiently in time domain. Its detailed description can be found in [8]. These methods were implemented in a number of computer codes, such as SPEAG SEMCAD™ (FDTD), CST Studio Suite™ (FIT, TLM) and many others.

The second group deals with Eqs. (14)–(17) for field phasors, supposing harmonic time dependence of fields. This group provides solution in Frequency (or Spectral) Domain (FD). The finite element method (FEM) is one of the most efficient representatives of this group. The first application of this method to the solution of mechanical problems refers to the year 1943 [9]. The book [10] contains detailed description of the method, which is rather universal and accurate. The Method of Moments (MoMs) and its varieties are also frequently used in FD. In contrast with FEM, MoM uses integral form of basic equations, where electric current density distribution on conducting surfaces excited by external sources is unknown. The book [11] reflects modern state of the MoM. Mentioned methods were implemented in codes ANSIS HFSS™ (FEM), Altair FEKO™ (MoM, FEM) and other commercial computer codes. Of course, there exist many other numerical methods and various implementations.

Delaunay tessellation is a common way of mesh building. It includes several stages. Firstly, a surface triangle mesh is generated. Then, a volume mesh based on the surface mesh and covering the whole computational region is being built. The method can build a rather good mesh without self-intersections. Some commercial codes, such as SYMMETRIX MeshSym™,

Computer Simulation of High‐Frequency Electromagnetic Fields

http://dx.doi.org/10.5772/67497

199

Unfortunately, this method cannot be applied to tessellation of complex geometrical models consisting of hundreds of parts made from various materials. Usually, an electromagnetic simulation code imports such models from various CAD systems (see, e.g. Figure 1). As a rule, imported models contain a number of errors caused by insufficient attention of a designer or arising in the process of graphic formats transforms. As a result, mesh generator fails to build the mesh.

On the other side, CAD models are often excessively detailed. They contain peculiarities, not influenced on electromagnetic field. Figure 2 shows an example of such extra detailed CAD model. The application of standard mesh generator to such a model can result in excess mesh size. Correction and simplification of the model manually needs several working days of the

An advanced method of mesh generation was developed and implemented in the RFS simulation code. The algorithm begins from assigning materials and attributes to the model parts. Most important objects, such as ports, printed circuit board (PCB) and antennas, are labelled as 'electrically important'. This is the only stage of the algorithm, which is made manually. This stage can be omitted for simple models. Then the code builds surface meshes for each detail of the system. The third stage presumes elimination of individual meshes interceptions and building united surface mesh. Electrically important, metal parts and parts with higher permittivity have priority in this process. At last, a global volume mesh is generated [13]. The user

can be used to generate mesh by this procedure.

qualified engineer.

Figure 1. A part of the handset CAD model.

Most of modern methods allow implementation both in TD and in FD.

Because it is impossible to give detailed representation of all numerical methods in a limited space, we give here more detailed information about the FEM method, which was implemented in the computer code radio frequency simulator (RFS) [12].

## 5. Finite element method for electromagnetics

#### 5.1. Main features of the method

The finite element method belongs to the variational methods of solving partial differential equation (PDE). It presumes formulating a functional, which is stationary (has minimum or maximum value) on the equation solution. In order to find functional extremum, the computational region is partitioned on a number of subregions (finite elements). After that, we approximate an unknown function in each finite element (FE) by superposition of basis functions. Basis functions have to be simple and form a linearly independent system. Applying Ritz method or Galerkin algorithm, we fulfil discretization of the problem, that is, transition from the partial difference equation to the system of linear algebraic equations (SLAE). Numerical solution of the SLAE gives unknown coefficients of basic functions, which are used to restore electromagnetic field. At last, needed components of electromagnetic field and system parameters are calculated.

#### 5.2. Mesh generation

Partitioning of the computational region (mesh generation) is the first stage of problem solution by FEM. It supposes dividing the computation region on a set of subregions—finite elements. FEs must densely fill the region and be nearly conformal to its border. In contrast to FDTD or FIT methods, a FEM mesh can be irregular and contain FEs of different forms. Tetrahedron FEs are used most commonly, because they allow dense packing and quite correctly approximate curvilinear borders. However, mesh generation for regions with complex forms filled with different materials is a very complex task.

Delaunay tessellation is a common way of mesh building. It includes several stages. Firstly, a surface triangle mesh is generated. Then, a volume mesh based on the surface mesh and covering the whole computational region is being built. The method can build a rather good mesh without self-intersections. Some commercial codes, such as SYMMETRIX MeshSym™, can be used to generate mesh by this procedure.

Unfortunately, this method cannot be applied to tessellation of complex geometrical models consisting of hundreds of parts made from various materials. Usually, an electromagnetic simulation code imports such models from various CAD systems (see, e.g. Figure 1). As a rule, imported models contain a number of errors caused by insufficient attention of a designer or arising in the process of graphic formats transforms. As a result, mesh generator fails to build the mesh.

On the other side, CAD models are often excessively detailed. They contain peculiarities, not influenced on electromagnetic field. Figure 2 shows an example of such extra detailed CAD model. The application of standard mesh generator to such a model can result in excess mesh size. Correction and simplification of the model manually needs several working days of the qualified engineer.

An advanced method of mesh generation was developed and implemented in the RFS simulation code. The algorithm begins from assigning materials and attributes to the model parts. Most important objects, such as ports, printed circuit board (PCB) and antennas, are labelled as 'electrically important'. This is the only stage of the algorithm, which is made manually. This stage can be omitted for simple models. Then the code builds surface meshes for each detail of the system. The third stage presumes elimination of individual meshes interceptions and building united surface mesh. Electrically important, metal parts and parts with higher permittivity have priority in this process. At last, a global volume mesh is generated [13]. The user

Figure 1. A part of the handset CAD model.

The second group deals with Eqs. (14)–(17) for field phasors, supposing harmonic time dependence of fields. This group provides solution in Frequency (or Spectral) Domain (FD). The finite element method (FEM) is one of the most efficient representatives of this group. The first application of this method to the solution of mechanical problems refers to the year 1943 [9]. The book [10] contains detailed description of the method, which is rather universal and accurate. The Method of Moments (MoMs) and its varieties are also frequently used in FD. In contrast with FEM, MoM uses integral form of basic equations, where electric current density distribution on conducting surfaces excited by external sources is unknown. The book [11] reflects modern state of the MoM. Mentioned methods were implemented in codes ANSIS HFSS™ (FEM), Altair FEKO™ (MoM, FEM) and other commercial computer codes. Of course, there exist many other numerical methods and various

Because it is impossible to give detailed representation of all numerical methods in a limited space, we give here more detailed information about the FEM method, which was

The finite element method belongs to the variational methods of solving partial differential equation (PDE). It presumes formulating a functional, which is stationary (has minimum or maximum value) on the equation solution. In order to find functional extremum, the computational region is partitioned on a number of subregions (finite elements). After that, we approximate an unknown function in each finite element (FE) by superposition of basis functions. Basis functions have to be simple and form a linearly independent system. Applying Ritz method or Galerkin algorithm, we fulfil discretization of the problem, that is, transition from the partial difference equation to the system of linear algebraic equations (SLAE). Numerical solution of the SLAE gives unknown coefficients of basic functions, which are used to restore electromagnetic field. At last, needed components of electromagnetic field and system param-

Partitioning of the computational region (mesh generation) is the first stage of problem solution by FEM. It supposes dividing the computation region on a set of subregions—finite elements. FEs must densely fill the region and be nearly conformal to its border. In contrast to FDTD or FIT methods, a FEM mesh can be irregular and contain FEs of different forms. Tetrahedron FEs are used most commonly, because they allow dense packing and quite correctly approximate curvilinear borders. However, mesh generation for regions with com-

Most of modern methods allow implementation both in TD and in FD.

implemented in the computer code radio frequency simulator (RFS) [12].

5. Finite element method for electromagnetics

plex forms filled with different materials is a very complex task.

5.1. Main features of the method

eters are calculated.

5.2. Mesh generation

implementations.

198 Computer Simulation

can control mesh quality (maximum to minimum ratio of tetrahedron's dimensions and other

We proved this algorithm on more than 190 models of mobile handsets and showed its 100% reliability. Generated meshes contained more than 1.5 billion tetrahedrons. Figure 3 shows a CAD model of the handset (a), surface mesh (b) and volume mesh (c), built by the described

Consider a region V, delimited by a surface S, where electromagnetic field has to be calculated. We suppose that the region V is filled by linear medium. Applying curl operator to Eq. (15) and

Ė ¼ �jkη0J

μ0=ε<sup>0</sup>

Computer Simulation of High‐Frequency Electromagnetic Fields

w<sup>m</sup> ¼ ζm1∇ζm<sup>2</sup> � ζm2∇ζm1, (20)

ζ<sup>p</sup> ¼ ap þ bpx þ cpy þ dpz, (21)

� imp, (18)

http://dx.doi.org/10.5772/67497

201

<sup>p</sup> <sup>¼</sup> <sup>120</sup><sup>π</sup> ohm is the intrinsic

∇ · Ė: (19)

�1=<sup>2</sup> is the light

<sup>r</sup> <sup>∇</sup> · <sup>Ė</sup>Þ þ <sup>j</sup>ση0k<sup>Ė</sup> � <sup>ε</sup>\_rk<sup>2</sup>

� imp is the imposed current density, <sup>k</sup> <sup>¼</sup> <sup>ω</sup>=<sup>c</sup> is wave number, <sup>c</sup> ¼ ðε0μ0<sup>Þ</sup>

H � <sup>¼</sup> <sup>j</sup> kη0μ\_ <sup>r</sup>

impedance of free space. After solving this equation, we can find magnetic field by means of

To solve Eq. (16), we divide the computation region on finite elements and approximate the field in each FE by superposition of basis functions. It is natural to approximate vector unknown function in Eq. (16) by a set of vector basis functions. We name such variant of FEM

Nedelec [14] proposed a set of vector basis functions for triangle and tetrahedral FE. These

where ζmn is the barycentric function of the vertex (node) node n of the edge m (m ¼ 1, …, 6, n ¼ 1, 2). The barycentric functions are polynomials of the first order. For a given

where p is vertex number. As a tetrahedron has six edges, Nedelec's basis consists of six functions. Figure 4 shows the numbering scheme for nodes and edges, needed for the correct

substituting into result Eq. (14), we get PDE of the second order

velocity in free space, <sup>σ</sup> is medium conductivity and <sup>η</sup><sup>0</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffi

functions are associated with tetrahedron edges and have the form

Coefficients in Eq. (17) can be found from the equation

<sup>∇</sup> · <sup>ð</sup>μ\_ �<sup>1</sup>

as vector finite element method (VFEM).

Descartes coordinate system

use of Nedelec's functions.

mesh parameters).

5.3. Basic equation

algorithm.

where J

Eq. (15):

Figure 2. A part of the CAD model with excessive details.

Figure 3. CAD model of the handset (a), united surface mesh (b) and volume mesh (c).

can control mesh quality (maximum to minimum ratio of tetrahedron's dimensions and other mesh parameters).

We proved this algorithm on more than 190 models of mobile handsets and showed its 100% reliability. Generated meshes contained more than 1.5 billion tetrahedrons. Figure 3 shows a CAD model of the handset (a), surface mesh (b) and volume mesh (c), built by the described algorithm.

#### 5.3. Basic equation

Figure 2. A part of the CAD model with excessive details.

200 Computer Simulation

Figure 3. CAD model of the handset (a), united surface mesh (b) and volume mesh (c).

Consider a region V, delimited by a surface S, where electromagnetic field has to be calculated. We suppose that the region V is filled by linear medium. Applying curl operator to Eq. (15) and substituting into result Eq. (14), we get PDE of the second order

$$\nabla \times (\dot{\mu}\_r^{-1} \nabla \times \dot{\mathbf{E}}) + j \sigma \eta\_0 k \dot{\mathbf{E}} - \dot{\varepsilon}\_r k^2 \dot{\mathbf{E}} = -jk \eta\_0 \dot{\mathbf{J}}^{imp},\tag{18}$$

where J � imp is the imposed current density, <sup>k</sup> <sup>¼</sup> <sup>ω</sup>=<sup>c</sup> is wave number, <sup>c</sup> ¼ ðε0μ0<sup>Þ</sup> �1=<sup>2</sup> is the light velocity in free space, <sup>σ</sup> is medium conductivity and <sup>η</sup><sup>0</sup> <sup>¼</sup> ffiffiffiffiffiffiffiffiffiffiffiffi μ0=ε<sup>0</sup> <sup>p</sup> <sup>¼</sup> <sup>120</sup><sup>π</sup> ohm is the intrinsic impedance of free space. After solving this equation, we can find magnetic field by means of Eq. (15):

$$
\dot{\mathbf{H}} = \frac{\dot{j}}{k\eta\_0 \dot{\mu}\_r} \nabla \times \dot{\mathbf{E}}.\tag{19}
$$

To solve Eq. (16), we divide the computation region on finite elements and approximate the field in each FE by superposition of basis functions. It is natural to approximate vector unknown function in Eq. (16) by a set of vector basis functions. We name such variant of FEM as vector finite element method (VFEM).

Nedelec [14] proposed a set of vector basis functions for triangle and tetrahedral FE. These functions are associated with tetrahedron edges and have the form

$$\mathbf{w}\_m = \zeta\_{m1}\nabla\zeta\_{m2} - \zeta\_{m2}\nabla\zeta\_{m1},\tag{20}$$

where ζmn is the barycentric function of the vertex (node) node n of the edge m (m ¼ 1, …, 6, n ¼ 1, 2). The barycentric functions are polynomials of the first order. For a given Descartes coordinate system

$$
\zeta\_p = a\_p + b\_p \mathbf{x} + c\_p \mathbf{y} + d\_p \mathbf{z},\tag{21}
$$

where p is vertex number. As a tetrahedron has six edges, Nedelec's basis consists of six functions. Figure 4 shows the numbering scheme for nodes and edges, needed for the correct use of Nedelec's functions.

Coefficients in Eq. (17) can be found from the equation

Figure 4. Numbering order of a tetrahedron's vertices and edges.

$$
\begin{bmatrix} b\_1 & c\_1 & d\_1 & a\_1 \\ b\_2 & c\_2 & d\_2 & a\_2 \\ b\_3 & c\_3 & d\_3 & a\_3 \\ b\_4 & c\_4 & d\_4 & a\_4 \end{bmatrix} = \begin{bmatrix} \mathbf{x}\_1 & \mathbf{x}\_2 & \mathbf{x}\_3 & \mathbf{x}\_4 \\ \mathbf{y}\_1 & \mathbf{y}\_2 & \mathbf{y}\_3 & \mathbf{y}\_4 \\ \mathbf{z}\_1 & \mathbf{z}\_2 & \mathbf{z}\_3 & \mathbf{z}\_4 \\ \mathbf{1} & \mathbf{1} & \mathbf{1} & \mathbf{1} \end{bmatrix}^{-1} \tag{22}
$$

4. Circulation of function w<sup>m</sup> along its edge is equal to unity.

exist basis functions of order 2.5 and higher, but they are rarely used.

<sup>ð</sup>rÞ ¼ <sup>X</sup> 6

> Lq Ė q

dron volume Vq. As a result, we obtain the so-called week form of Eq. (16):

<sup>m</sup>ÞdV <sup>þ</sup> jk<sup>ð</sup>

<sup>m</sup>�d<sup>S</sup> � jk<sup>ð</sup>

<sup>σ</sup>E<sup>q</sup> � <sup>w</sup><sup>q</sup>

Vq <sup>η</sup>0σw<sup>q</sup>

Vq η0J � imp � <sup>w</sup><sup>q</sup>

ð Vq

n¼1 x\_ q nw<sup>q</sup>

along edges and in normal direction to them.

intensity in q-th tetrahedron is approximated as

where Nq is the total number of tetrahedrons.

due to the fourth property of basis functions.

Þ� � <sup>w</sup><sup>q</sup>

<sup>m</sup>dV þ jη0k

approximation (19) gives a system of linear algebraic equations:

<sup>n</sup>Þ�ð<sup>∇</sup> · <sup>w</sup><sup>q</sup>

Vq <sup>m</sup> ¼ �<sup>ð</sup>

Coefficients x\_ <sup>q</sup>

along an edge

basis function w<sup>q</sup>

<sup>½</sup><sup>∇</sup> · <sup>ð</sup>μ\_ �<sup>1</sup>

¼ �∮ Sq

½ðμ\_ �<sup>1</sup> <sup>r</sup> ∇ · Ė q <sup>Þ</sup> · <sup>w</sup><sup>q</sup>

<sup>r</sup> <sup>∇</sup> · <sup>E</sup><sup>q</sup>

ð Vq

> X 6

n¼1 x\_ q n½ ð Vq <sup>ð</sup>μ\_ �<sup>1</sup> <sup>r</sup> <sup>∇</sup> · <sup>w</sup><sup>q</sup> Ė q

Full first-order polynomial set for electric field approximation contains 12 functions (3 E-vector projections, each needs four functions). However, Nedelec's set contains only six basis functions. Hence, it does not form full first-order polynomial system. The order of approximation by these functions is lower than one. Finite elements with these functions are called FE of order 1/2 or CT/LN (constant tangential/linear normal) according to the behaviour of these functions

Striving for increasing approximation accuracy led to creating another set of basis functions with higher order of approximation. A set of 20 functions (12 associated with edges and eight with faces) forms FE of order 1.5 or LN/QN (linear tangential/quadratic normal) [15, 16]. There

With the aim of simplicity, we consider further only Nedelec's basis functions. So, electric field

� <sup>d</sup>l<sup>m</sup> ¼ �<sup>X</sup>

Let us consider q-th tetrahedron. According to Galerkin's method, we multiply Eq. (16) by

<sup>m</sup>dV � k 2 ð Vq

<sup>n</sup> � <sup>w</sup><sup>q</sup>

<sup>m</sup>dV � <sup>k</sup><sup>2</sup>

Transformation of the first term of this equation and subsequent substitution into result

6

n¼1 x\_ q nw<sup>q</sup>

<sup>m</sup> have the sense of voltage on edge (with reverse sign). Really, taking integral

m, belonging to this tetrahedron and integrate the product over the tetrahe-

<sup>n</sup>dl<sup>m</sup> ¼ �x\_ <sup>q</sup>

<sup>ε</sup>\_ <sup>r</sup>E<sup>q</sup> � <sup>w</sup><sup>q</sup>

ð Vq ε\_ <sup>r</sup>w<sup>q</sup> <sup>n</sup> � <sup>w</sup><sup>q</sup> <sup>m</sup>dV�

<sup>m</sup>dV ¼ �jkη<sup>0</sup>

<sup>m</sup>dV, m ¼ 1, 2, …, 6; q ¼ 1, 2, …, Nq,

<sup>n</sup>ðrÞ, q ¼ 1, …, Nq, (23)

Computer Simulation of High‐Frequency Electromagnetic Fields

http://dx.doi.org/10.5772/67497

203

<sup>m</sup>, (24)

ð Vq J imp � <sup>w</sup><sup>q</sup>

<sup>m</sup>dV:

(25)

where xp, yp, zp, p ¼ 1…4 are the coordinates of the p-th tetrahedron's node.

Nedelec's basis functions have the following properties:


4. Circulation of function w<sup>m</sup> along its edge is equal to unity.

Full first-order polynomial set for electric field approximation contains 12 functions (3 E-vector projections, each needs four functions). However, Nedelec's set contains only six basis functions. Hence, it does not form full first-order polynomial system. The order of approximation by these functions is lower than one. Finite elements with these functions are called FE of order 1/2 or CT/LN (constant tangential/linear normal) according to the behaviour of these functions along edges and in normal direction to them.

Striving for increasing approximation accuracy led to creating another set of basis functions with higher order of approximation. A set of 20 functions (12 associated with edges and eight with faces) forms FE of order 1.5 or LN/QN (linear tangential/quadratic normal) [15, 16]. There exist basis functions of order 2.5 and higher, but they are rarely used.

With the aim of simplicity, we consider further only Nedelec's basis functions. So, electric field intensity in q-th tetrahedron is approximated as

$$\dot{\mathbf{E}}^{\mathcal{q}}(\mathbf{r}) = \sum\_{n=1}^{6} \dot{\mathbf{x}}\_{n}^{\mathcal{q}} \mathbf{w}\_{n}^{\mathcal{q}}(\mathbf{r}), \quad \boldsymbol{\eta} = 1, \ldots, N\_{\boldsymbol{\eta}}, \tag{23}$$

where Nq is the total number of tetrahedrons.

Coefficients x\_ <sup>q</sup> <sup>m</sup> have the sense of voltage on edge (with reverse sign). Really, taking integral along an edge

$$\mathcal{V}\_{m}^{q} = -\int\_{L\_{q}} \dot{\mathbf{E}}^{q} \cdot d\mathbf{l}\_{m} = -\sum\_{n=1}^{6} \dot{\mathbf{x}}\_{n}^{q} \mathbf{w}\_{n}^{q} d\mathbf{l}\_{m} = -\dot{\mathbf{x}}\_{m}^{q},\tag{24}$$

due to the fourth property of basis functions.

b<sup>1</sup> c<sup>1</sup> d<sup>1</sup> a<sup>1</sup> b<sup>2</sup> c<sup>2</sup> d<sup>2</sup> a<sup>2</sup> b<sup>3</sup> c<sup>3</sup> d<sup>3</sup> a<sup>3</sup> b<sup>4</sup> c<sup>4</sup> d<sup>4</sup> a<sup>4</sup>

where xp, yp, zp, p ¼ 1…4 are the coordinates of the p-th tetrahedron's node.

.

2. They have zero divergence and their curl is constant inside the FE. This property provides

3. Projection of the function w<sup>m</sup> on the edge m is constant, and its projection on the other edges is equal to zero. Hence, these functions provide continuity of tangential electric field

x<sup>1</sup> x<sup>2</sup> x<sup>3</sup> x<sup>4</sup> y<sup>1</sup> y<sup>2</sup> y<sup>3</sup> y<sup>4</sup> z<sup>1</sup> z<sup>2</sup> z<sup>3</sup> z<sup>4</sup> 1111

�1

: (22)

Figure 4. Numbering order of a tetrahedron's vertices and edges.

202 Computer Simulation

Nedelec's basis functions have the following properties:

the absence of spurious solutions of the problem.

on the border between neighbouring tetrahedrons.

1. Dimension of these functions is L�<sup>1</sup>

Let us consider q-th tetrahedron. According to Galerkin's method, we multiply Eq. (16) by basis function w<sup>q</sup> m, belonging to this tetrahedron and integrate the product over the tetrahedron volume Vq. As a result, we obtain the so-called week form of Eq. (16):

$$\int\_{V\_q} \left[\nabla \times (\dot{\mu}\_r^{-1} \nabla \times \mathbf{E}^q) \right] \cdot \mathbf{w}\_m^q dV + j\eta\_0 k \int\_{V\_q} \sigma \mathbf{E}^q \cdot \mathbf{w}\_m^q dV - k^2 \int\_{V\_q} \dot{\varepsilon}\_r \mathbf{E}^q \cdot \mathbf{w}\_m^q dV = -jk\eta\_0 \int\_{V\_q} \mathbf{J}^{imp} \cdot \mathbf{w}\_m^q dV.$$

Transformation of the first term of this equation and subsequent substitution into result approximation (19) gives a system of linear algebraic equations:

$$\begin{split} \sum\_{n=1}^{6} \dot{\mathbf{x}}\_{n}^{q} [\int\_{V\_{q}} (\dot{\mu}\_{r}^{-1} \nabla \times \mathbf{w}\_{n}^{q}) \cdot (\nabla \times \mathbf{w}\_{m}^{q}) dV + j \mathbf{k} \int\_{V\_{q}} \eta\_{0} \sigma \mathbf{w}\_{n}^{q} \cdot \mathbf{w}\_{m}^{q} dV - \boldsymbol{k}^{2} \int\_{V\_{q}} \dot{\varepsilon}\_{r} \mathbf{w}\_{n}^{q} \cdot \mathbf{w}\_{m}^{q} dV] \\ = - \oint\_{S\_{q}} [(\dot{\mu}\_{r}^{-1} \nabla \times \dot{\mathbf{E}}^{\boldsymbol{\theta}}) \times \mathbf{w}\_{m}^{\boldsymbol{\theta}}] d\mathbf{S} - j \boldsymbol{k} \int\_{V\_{q}} \eta\_{0} \dot{\mathbf{J}}^{\text{imp}} \cdot \mathbf{w}\_{m}^{q} dV, \quad m = 1, 2, \dots, 6; \quad q = 1, 2, \dots, N\_{\boldsymbol{\theta}}, \tag{25} \end{split} \tag{26}$$

where Nq is the total number of FE in the computation region.

The system of Eq. (25) can be written in matrix form:

$$\overline{\mathbf{Q}}^{\dagger}\mathbf{X}^{q} \equiv (\overline{\mathbf{T}}^{\dagger} + jk\overline{\mathbf{U}}^{\dagger} - k^{2}\overline{\mathbf{R}}^{q})\overline{\mathbf{X}}^{q} = -\overline{\mathbf{S}}^{q} - jk\overline{\mathbf{B}}^{q},\tag{26}$$

• We introduce matrix G with components

• Integrals (32) are calculated:

dimensions instead of 12.

numbers from 1 to Nc ¼ 6: 2. Choose the tetrahedron with q = 2

Figure 5. Two neighbour tetrahedrons. Edges e11 and e22

v4<sup>1</sup> and v4<sup>2</sup> are common for both FEs.

In these expressions, V is the volume of the tetrahedron.

The code RFS uses the next algorithm for building the global matrix:

gij <sup>¼</sup> <sup>V</sup>�<sup>1</sup>

<sup>G</sup> <sup>¼</sup> <sup>1</sup> 20

Using integration formula for three-dimensional (3D) simplex, we can prove that

� � � � � � � � ð V

2111 1211 1121 1112

hmn ¼ Vðϕm2,n2gm1,n<sup>1</sup> � ϕm2,n1gm1,n<sup>2</sup> � ϕm1,n2gm2,n<sup>1</sup> þ ϕm1,n1gm2,n2Þ, m, n ¼ 1, …, 6: (38)

Construction of the global matrix (assembling) for all tetrahedrons is a rather complex task, because many nodes, edges and faces are common for neighbour FEs (see Figure 5, which shows two neighbour tetrahedrons, separated for reader's convenience). They have three common vertices and three common edges. Hence, the united matrix for two FEs has only 9

1. Calculate a local matrix for the tetrahedron number one (q = 1). Give its edges global

, e31 and e32

, e5<sup>1</sup> and e6<sup>2</sup> and vertices v1<sup>1</sup> and v12

, v2<sup>1</sup> and v3<sup>2</sup>

,

� � � � � � � �

f mn ¼ 4VðVm1,m<sup>2</sup> � Vn1,n2ÞdV; (37)

ζiζjdV: (35)

Computer Simulation of High‐Frequency Electromagnetic Fields

: (36)

http://dx.doi.org/10.5772/67497

205

where Q<sup>q</sup> , T<sup>q</sup> , U<sup>q</sup> , R<sup>q</sup> are square matrices having dimension 6 · 6, S q , B<sup>q</sup> are column vectors having six components, dependent on the given boundary conditions and implicit current density in the q-th finite element. Depending on BC, S<sup>q</sup> can sometimes be a square 6 · 6 matrix.

Equation (26) is called the local matrix equation, and matrix Q<sup>q</sup> is the local matrix of q-th FE. Such matrices should be built for every FE in the computational region.

Components of matrices T<sup>q</sup> , U<sup>q</sup> , R<sup>q</sup> and vectors S q , B<sup>q</sup> are defined by the expressions

$$t\_{mn} = \dot{\mu}\_r^{-1} \int\_V (\nabla \times \mathbf{w}\_m) \cdot (\nabla \times \mathbf{w}\_n) dV;\tag{27}$$

$$
\mu\_{\rm mn} = \eta\_0 \sigma \int\_V \mathbf{w}\_{\rm m} \cdot \mathbf{w}\_{\rm n} dV;\tag{28}
$$

$$r\_{mn} = \dot{\varepsilon}\_r \int\_V \mathbf{w}\_n \cdot \mathbf{w}\_m dV;\tag{29}$$

$$s\_m = \dot{\mu}\_r^{-1} \oint\_S [(\nabla \times \dot{\mathbf{E}}) \times \mathbf{w}\_m] d\mathbf{S};\tag{30}$$

$$b\_m = \eta\_0 \int\_V \mathbf{j}^{imp} \cdot \mathbf{w}\_m dV. \tag{31}$$

Medium in these expressions is supposed to be homogeneous inside the FE. Index 'q' denoted FE number is omitted. We can see that in order to calculate elements (27)–(31), it is sufficient to calculate two integrals:

$$f\_{mn} = \int\_{V} (\nabla \times \mathbf{w}\_n) \cdot (\nabla \times \mathbf{w}\_m) dV \text{ and } h\_{mn} = \int\_{V} (\mathbf{w}\_n \cdot \mathbf{w}\_m) dV. \tag{32}$$

We use the next algorithm for the calculation of matrices and vectors elements:

#### For each tetrahedron


$$\mathbf{v}\_{i\uparrow} = \nabla \zeta\_i \times \nabla \zeta\_{\uparrow} = (c\_i d\_{\uparrow} - c\_{\rangle} d\_i)\widehat{\mathbf{x}} + (d\_i b\_{\slash} - d\_{\rangle} b\_i)\widehat{\mathbf{y}} + (b\_i c\_{\slash} - b\_{\rangle} c\_i)\widehat{\mathbf{z}};\tag{33}$$

$$
\varphi\_{ij} = \nabla \zeta\_i \cdot \nabla \zeta\_j = b\_i b\_j + c\_i c\_j + d\_i d\_j. \tag{34}
$$

It can be easily shown that ∇ · w<sup>m</sup> ¼ 2vm1,m2.

• We introduce matrix G with components

where Nq is the total number of FE in the computation region.

<sup>X</sup><sup>q</sup> � ðT<sup>q</sup> <sup>þ</sup> jkU<sup>q</sup> � <sup>k</sup>

Such matrices should be built for every FE in the computational region.

tmn <sup>¼</sup> <sup>μ</sup>\_ �<sup>1</sup> r ð V

, R<sup>q</sup> and vectors S

umn ¼ η0σ

rmn ¼ ε\_<sup>r</sup>

bm ¼ η<sup>0</sup>

sm <sup>¼</sup> <sup>μ</sup>\_ �<sup>1</sup>

, U<sup>q</sup>

2 Rq

having six components, dependent on the given boundary conditions and implicit current density in the q-th finite element. Depending on BC, S<sup>q</sup> can sometimes be a square 6 · 6

Equation (26) is called the local matrix equation, and matrix Q<sup>q</sup> is the local matrix of q-th FE.

ð V

ð V

ð V J � imp

ð∇ · wnÞ�ð∇ · wmÞdV and hmn ¼

We use the next algorithm for the calculation of matrices and vectors elements:

• Coefficients of the barycentric functions are calculated by means of Eq. (22).

vij ¼ ∇ζ<sup>i</sup> · ∇ζ<sup>j</sup> ¼ ðcidj � cjdiÞx

It can be easily shown that ∇ · w<sup>m</sup> ¼ 2vm1,m2.

Medium in these expressions is supposed to be homogeneous inside the FE. Index 'q' denoted FE number is omitted. We can see that in order to calculate elements (27)–(31), it is sufficient to

q

, R<sup>q</sup> are square matrices having dimension 6 · 6, S

<sup>Þ</sup>X<sup>q</sup> ¼ �<sup>S</sup>

q � jkB<sup>q</sup>

q

, B<sup>q</sup> are defined by the expressions

ð∇ · wmÞ�ð∇ · wnÞdV; (27)

w<sup>m</sup> � wndV; (28)

w<sup>n</sup> � wmdV; (29)

� wmdV: (31)

ðw<sup>n</sup> � wmÞdV: (32)

\_; (33)

<sup>r</sup> ∮ <sup>S</sup>½ð∇ · ĖÞ · wm�dS; (30)

\_ þ ðbicj � bjciÞ<sup>z</sup>

ϕij ¼ ∇ζ<sup>i</sup> � ∇ζ<sup>j</sup> ¼ bibj þ cicj þ didj: (34)

ð V

\_ þ ðdibj � djbiÞ<sup>y</sup>

, (26)

, B<sup>q</sup> are column vectors

The system of Eq. (25) can be written in matrix form: Qq

where Q<sup>q</sup>

204 Computer Simulation

matrix.

, T<sup>q</sup> , U<sup>q</sup>

Components of matrices T<sup>q</sup>

calculate two integrals:

For each tetrahedron

f mn ¼ ð V

• Auxiliary values are calculated:

$$\mathcal{g}\_{i\dot{\jmath}} = V^{-1} \int\_{V} \zeta\_i \zeta\_{\dot{\jmath}} dV. \tag{35}$$

Using integration formula for three-dimensional (3D) simplex, we can prove that

$$\overline{\mathbf{G}} = \frac{1}{20} \begin{vmatrix} 2 & 1 & 1 & 1 \\ 1 & 2 & 1 & 1 \\ 1 & 1 & 2 & 1 \\ 1 & 1 & 1 & 2 \end{vmatrix} . \tag{36}$$

• Integrals (32) are calculated:

$$f\_{mn} = 4V(\mathbf{V}\_{m1,m2} \cdot \mathbf{V}\_{n1,n2})dV;\tag{37}$$

$$h\_{\rm nn} = V(\boldsymbol{\varphi}\_{\rm m2,n2}\mathbf{g}\_{\rm m1,n1} - \boldsymbol{\varphi}\_{\rm m2,n1}\mathbf{g}\_{\rm m1,n2} - \boldsymbol{\varphi}\_{\rm m1,n2}\mathbf{g}\_{\rm m2,n1} + \boldsymbol{\varphi}\_{\rm m1,n1}\mathbf{g}\_{\rm m2,n2}), \quad m, n = 1, \ldots, 6. \tag{38}$$

In these expressions, V is the volume of the tetrahedron.

Construction of the global matrix (assembling) for all tetrahedrons is a rather complex task, because many nodes, edges and faces are common for neighbour FEs (see Figure 5, which shows two neighbour tetrahedrons, separated for reader's convenience). They have three common vertices and three common edges. Hence, the united matrix for two FEs has only 9 dimensions instead of 12.

The code RFS uses the next algorithm for building the global matrix:


Figure 5. Two neighbour tetrahedrons. Edges e11 and e22 , e31 and e32 , e5<sup>1</sup> and e6<sup>2</sup> and vertices v1<sup>1</sup> and v12 , v2<sup>1</sup> and v3<sup>2</sup> , v4<sup>1</sup> and v4<sup>2</sup> are common for both FEs.

a. Choose the edge of this tetrahedron with m = 1. Check whether this edge coincides with any edge of the previous tetrahedron. If it does, calculate its matrix elements and add them to corresponding matrix elements of the previous tetrahedron. In the other case, give the edge global number Nc þ 1 and fill a corresponding element of the global matrix. Set Nc ¼ Nc þ 1.

• Impedance surface. Impedance surface has finite surface conductivity σs. Such surface is a good approximation of the metal surface in microwave frequency band. Using Leontovich

> ð S

These values must be added to the corresponding elements of the local matrix. Of course, they are nonzero only for functions associated with edges lying on the impedance surface.

Together with wave ports, lumped ports (LPs) are often used. These ports are defined by current I and intrinsic admittance Zi. LPs are used when excitation source dimensions are very low compared to the whole-system dimensions or wavelength. Geometrical model of LP is a segment of a straight line (linear LP) or a rectangle (planar LP). Linear and planar LPs are implemented in the RFS code. Arbitrary number of tetrahedron's edges or (and) faces can cover such LPs. We suppose that surface current density and electric field distributions on the PL are homogeneous. Values of right-hand vector for a planar LP are

model of the LP excludes limitations imposed on the mesh generator by other models. • Lumped elements. Embedding of lumped elements (LEs) into 3D field model enlarges the capability of the code, as an LE excludes the necessity of generating very fine mesh inside such parts as resistors, capacitors and inductors. In the RFS code, linear and planar LEs

The presence of LEs modifies the global matrix by adding to diagonal matrix elements

where Y is the admittance of the LE, l is its length. The RFS code implements models of a resistor, capacitor and inductor, and their parallel connection. The number of LEs in the

The finite element method in frequency domain solves the problem at one specific frequency. Calculation of the system frequency response (e.g. S-matrix) needs solving the problem many times (often several hundreds or even thousands). This procedure takes much computation time.

• Excitation sources. Electromagnetic field in the system excites by electric current or by electric or magnetic fields, defined on a part of the system border. In the first case, we use Eq. (26) to calculate right-hand vector elements. In the second case, we define the so-called ports—surfaces on which nonzero excitation field exists. Frequently, we define port as cross section of the regular transmission line (TL). Such ports are called wave ports. If TL cross section has simple form (rectangular, circular, etc.), we can define field on the port analytically. Otherwise, we have to solve two-dimensional (2D) boundary problem.

w<sup>m</sup> � wndS: (41)

Computer Simulation of High‐Frequency Electromagnetic Fields

http://dx.doi.org/10.5772/67497

207

¼ I=w, I is port current and w is port width. Implemented

mm ¼ �jk0η0Yl, (42)

smn ¼ �jk0η0σ<sup>s</sup>

BC (10), we can derive expression

calculated by Eq. (31), where J

are implemented [17].

model is unlimited.

5.5. Fast-frequency sweep

terms

� imp

q0


Similar assembling procedure is valid for basis functions of higher order.

As a result, we get a global matrix of the problem, a global vector of the right side of Eq. (26) and a global vector of unknowns. Dimension of the global matrix is equal to the total number of coefficients (degrees of freedom) in formulas (23).

#### 5.4. Numerical approximation of boundary conditions


$$\mathbf{s}\_{m} = j k\_{0} \eta\_{0} \int\_{S} (\mathbf{H} \times \mathbf{w}\_{m}) \mathbf{e}\_{n} = 0. \tag{39}$$

This expression is equal to zero as the mixed product under the integral has two vectors (H, en) with the same direction. Therefore, PMC BC does not change local matrix.

• Absorption boundary conditions. Formula for matrix elements smn can be derived from Eq. (11):

$$\mathbf{s}\_{mn} = \mathbf{j}\mathbf{k}\_0 \sqrt{\mu\_r/\varepsilon\_r} \int\_S (\mathbf{e}\_n \times \mathbf{w}\_m) \cdot (\mathbf{e}\_n \times \mathbf{w}\_n) d\mathbf{S}.\tag{40}$$

Expressions for the higher order ABC can be found in [10].

• Impedance surface. Impedance surface has finite surface conductivity σs. Such surface is a good approximation of the metal surface in microwave frequency band. Using Leontovich BC (10), we can derive expression

$$\mathbf{s}\_{mn} = -jk\_0 \eta\_0 \sigma\_s \int\_{\mathcal{S}} \mathbf{w}\_m \cdot \mathbf{w}\_m dS. \tag{41}$$

These values must be added to the corresponding elements of the local matrix. Of course, they are nonzero only for functions associated with edges lying on the impedance surface.

• Excitation sources. Electromagnetic field in the system excites by electric current or by electric or magnetic fields, defined on a part of the system border. In the first case, we use Eq. (26) to calculate right-hand vector elements. In the second case, we define the so-called ports—surfaces on which nonzero excitation field exists. Frequently, we define port as cross section of the regular transmission line (TL). Such ports are called wave ports. If TL cross section has simple form (rectangular, circular, etc.), we can define field on the port analytically. Otherwise, we have to solve two-dimensional (2D) boundary problem.

Together with wave ports, lumped ports (LPs) are often used. These ports are defined by current I and intrinsic admittance Zi. LPs are used when excitation source dimensions are very low compared to the whole-system dimensions or wavelength. Geometrical model of LP is a segment of a straight line (linear LP) or a rectangle (planar LP). Linear and planar LPs are implemented in the RFS code. Arbitrary number of tetrahedron's edges or (and) faces can cover such LPs. We suppose that surface current density and electric field distributions on the PL are homogeneous. Values of right-hand vector for a planar LP are calculated by Eq. (31), where J � imp ¼ I=w, I is port current and w is port width. Implemented model of the LP excludes limitations imposed on the mesh generator by other models.

• Lumped elements. Embedding of lumped elements (LEs) into 3D field model enlarges the capability of the code, as an LE excludes the necessity of generating very fine mesh inside such parts as resistors, capacitors and inductors. In the RFS code, linear and planar LEs are implemented [17].

The presence of LEs modifies the global matrix by adding to diagonal matrix elements terms

$$\mathbf{q'}\_{mm} = -j\mathbf{k}\_0 \boldsymbol{\eta}\_0 \mathbf{Y} \mathbf{l},\tag{42}$$

where Y is the admittance of the LE, l is its length. The RFS code implements models of a resistor, capacitor and inductor, and their parallel connection. The number of LEs in the model is unlimited.

#### 5.5. Fast-frequency sweep

a. Choose the edge of this tetrahedron with m = 1. Check whether this edge coincides with any edge of the previous tetrahedron. If it does, calculate its matrix elements and add them to corresponding matrix elements of the previous tetrahedron. In the other case, give the edge global number Nc þ 1 and fill a corresponding element of the

3. Repeat point 2 for every tetrahedron in the region testing whether a given edge of the

As a result, we get a global matrix of the problem, a global vector of the right side of Eq. (26) and a global vector of unknowns. Dimension of the global matrix is equal to the total number

• Border between two dielectrics. A finite element mesh is built so that every tetrahedron lies inside one of the media. Hence, the common face of two neighbour tetrahedrons lies on the separation surface or approximate it in the case of curvilinear border. Due to basis function properties, tangential electric field of both tetrahedrons is equal on the face edges, consequently electric fields are equal on the whole face area. As a result, BCs (7) for electric field are satisfied automatically. As for magnetic field intensity, its continuity

• Electrical wall. According to Eq. (8), tangential electric field on PEC surface is equal to zero. Hence, three coefficients of Eq. (25) belonging to edges, lying on the surface, are equal to zero. Corresponding rows and columns of local matrix of the given tetrahedron

• Magnetic wall. According to Eq. (9), surface tangential magnetic field on the PMC is equal

This expression is equal to zero as the mixed product under the integral has two vectors

ðH · wmÞe<sup>n</sup> ¼ 0: (39)

ðe<sup>n</sup> · wmÞ�ðe<sup>n</sup> · wnÞdS: (40)

on the border between media with different permeabilities is not guaranteed.

to zero. We can see from Eqs. (27)–(31) that only term sm has to be changed:

ð S

(H, en) with the same direction. Therefore, PMC BC does not change local matrix.

S

ffiffiffiffiffiffiffiffiffiffiffi μr=ε<sup>r</sup> q ð

• Absorption boundary conditions. Formula for matrix elements smn can be derived from

sm ¼ jk0η<sup>0</sup>

smn ¼ jk<sup>0</sup>

Expressions for the higher order ABC can be found in [10].

tetrahedron coincides with any edge, previously included in analysis.

Similar assembling procedure is valid for basis functions of higher order.

global matrix. Set Nc ¼ Nc þ 1.

206 Computer Simulation

of coefficients (degrees of freedom) in formulas (23).

have to be eliminated.

Eq. (11):

5.4. Numerical approximation of boundary conditions

b. Do this procedure for all edges of the tetrahedron.

The finite element method in frequency domain solves the problem at one specific frequency. Calculation of the system frequency response (e.g. S-matrix) needs solving the problem many times (often several hundreds or even thousands). This procedure takes much computation time. Considerable economy of computational resources can be achieved by implementing the socalled fast-frequency sweep (FFS). The algorithm of FFS is based on model-order reduction (MOR) technique. According to the MOR, we seek solution in one frequency point and then expand left and right sides of Eq. (21) in frequency power series. Terms with the same power are equated and full solution as a function of frequency is restored. The RFS code implements one of such algorithms—well conditioned asymptotic wave form evaluation (WCAWE) [18]. This method was adapted to the RFS procedures, particularly for LE models, and showed high efficiency, making possible to calculate system fields and parameters in three-octave frequency band using full solution only in one frequency point. Hence, computation time for wideband systems decreases by tens and hundreds times without losing accuracy.

Computation results together with analytical data are shown in Table 1. The table also contains relative calculation error δ. Azimuthal inhomogeneous modes are twice degenerated. Numerical errors remove this degeneration, so two eigen frequencies for each mode are given. Degenerated modes differ by their field, turned by 90° in azimuth direction and computation results confirm this phenomenon. Difference between degenerated eigen frequencies is no

6768 15,026 26,580

13.3132 13.3139

17.6853 17.6877

18.3081 18.3102

EF, GHz δ, % EF, GHz δ, % EF, GHz δ, %

0.084 0.089

Computer Simulation of High‐Frequency Electromagnetic Fields

0.086 0.102

0.131 0.143 13.3104 13.311

http://dx.doi.org/10.5772/67497

17.6802 17.6811

18.299 18.291 0.063 0.067 209

0.058 0.063

0.082 0.083

Mesh size

13.319

17.705

18.3287

Table 1. Comparison of analytical and computational eigen frequencies for a cylindrical cavity.

The table demonstrates that eigen frequency calculation error decreases with mesh refinement

TM010 11.474 11.502 0.235 11.4905 0.135 11.4848 0.085

TM011 15.216 15.232 0.098 15.2287 0.077 15.2241 0.052

0.147 0.198

0.222 0.244

0.11 0.135

Another example, an electric dipole, consisted of two perfectly conducting cylinders, connected to the linear-lumped source (Figure 6). To solve this outer problem, we ought to constrain

Figure 6. A dipole model: a, dipole structure; b, mesh; c, 3D directivity chart; d, 2D directivity chart in the plane ϕ ¼ 0.

more than hundreds of per cent.

frequency (EF), GHz

TE111 13.302 13.317

TE211 17.670 17.696

TM110 18.284 18.3246

Mode Analytical eigen

and not exceeds 0.1% on the last mesh.

## 6. The RFS code description

The RFS code is designated to solve complex 3D electromagnetic problems, especially mobile handsets modelling. It contains graphic user interface for creating and editing geometric objects or for import CAD models. Geometric primitives, besides boxes, cylinders, spheres, pyramids and cones, include coaxial and strip lines and wizards for creating human head and hand models (phantoms), spiral antennas and other objects. All Boolean operations, such as extrusion, skinning and others, are also implemented. The code includes a rich material library.

Two types of meshing are available: exact automatic meshing based on Delaunay tessellation and hand meshing with automatic geometry corrections. We recommend the last one for meshing complex models, imported from CAD codes, for example, for full models of a handset. Two types of basis functions can be used: low order (CT/LN) and high order (LT/QN). The last is used as default, but for complex models with many small details, we recommend low-order basis.

The code solves several types of problems: driven solution in one frequency point, driven solution in frequency band with given frequency step (solution in each point or FFS), eigenmode, EMS and EMI problems. The code can perform parametric solution when one or several geometric or (and)material parameters changes in a given order. The results can be represented as functions of these parameters.

Post-processing includes the calculation of S-, Y- and Z-parameters, field distribution in a volume, on a surface or along a line, the calculation of cavities intrinsic impedance, specific absorption rate (SAR) and EMC/EMI parameters. Results are represented in one-dimensional (1D), 2D and 3D graphs, tables and can be exported in a file.

## 7. Computation results

Firstly, we validated the code by solving problems having analytical solution. One of such problems is an eigenvalue problem for cylindrical cavity. A cylindrical cavity has curvilinear surface, so calculation can demonstrate the quality of surface approximation by a tetrahedron mesh. The analysed cavity had radius 10 mm and height 15 mm. Cavity walls were supposed to be perfectly conducting. Eight eigen frequencies and eigen fields were calculated using high-order basis functions (LT/QN).

Computation results together with analytical data are shown in Table 1. The table also contains relative calculation error δ. Azimuthal inhomogeneous modes are twice degenerated. Numerical errors remove this degeneration, so two eigen frequencies for each mode are given. Degenerated modes differ by their field, turned by 90° in azimuth direction and computation results confirm this phenomenon. Difference between degenerated eigen frequencies is no more than hundreds of per cent.

Considerable economy of computational resources can be achieved by implementing the socalled fast-frequency sweep (FFS). The algorithm of FFS is based on model-order reduction (MOR) technique. According to the MOR, we seek solution in one frequency point and then expand left and right sides of Eq. (21) in frequency power series. Terms with the same power are equated and full solution as a function of frequency is restored. The RFS code implements one of such algorithms—well conditioned asymptotic wave form evaluation (WCAWE) [18]. This method was adapted to the RFS procedures, particularly for LE models, and showed high efficiency, making possible to calculate system fields and parameters in three-octave frequency band using full solution only in one frequency point. Hence, computation time for wideband

The RFS code is designated to solve complex 3D electromagnetic problems, especially mobile handsets modelling. It contains graphic user interface for creating and editing geometric objects or for import CAD models. Geometric primitives, besides boxes, cylinders, spheres, pyramids and cones, include coaxial and strip lines and wizards for creating human head and hand models (phantoms), spiral antennas and other objects. All Boolean operations, such as extrusion, skinning and others, are also implemented. The code includes a rich material library. Two types of meshing are available: exact automatic meshing based on Delaunay tessellation and hand meshing with automatic geometry corrections. We recommend the last one for meshing complex models, imported from CAD codes, for example, for full models of a handset. Two types of basis functions can be used: low order (CT/LN) and high order (LT/QN). The last is used as default, but for complex models with many small details, we recommend low-order basis.

The code solves several types of problems: driven solution in one frequency point, driven solution in frequency band with given frequency step (solution in each point or FFS), eigenmode, EMS and EMI problems. The code can perform parametric solution when one or several geometric or (and)material parameters changes in a given order. The results can be represented as functions of these parameters. Post-processing includes the calculation of S-, Y- and Z-parameters, field distribution in a volume, on a surface or along a line, the calculation of cavities intrinsic impedance, specific absorption rate (SAR) and EMC/EMI parameters. Results are represented in one-dimensional (1D), 2D

Firstly, we validated the code by solving problems having analytical solution. One of such problems is an eigenvalue problem for cylindrical cavity. A cylindrical cavity has curvilinear surface, so calculation can demonstrate the quality of surface approximation by a tetrahedron mesh. The analysed cavity had radius 10 mm and height 15 mm. Cavity walls were supposed to be perfectly conducting. Eight eigen frequencies and eigen fields were calculated using high-order basis func-

systems decreases by tens and hundreds times without losing accuracy.

6. The RFS code description

208 Computer Simulation

and 3D graphs, tables and can be exported in a file.

7. Computation results

tions (LT/QN).


Table 1. Comparison of analytical and computational eigen frequencies for a cylindrical cavity.

The table demonstrates that eigen frequency calculation error decreases with mesh refinement and not exceeds 0.1% on the last mesh.

Another example, an electric dipole, consisted of two perfectly conducting cylinders, connected to the linear-lumped source (Figure 6). To solve this outer problem, we ought to constrain

Figure 6. A dipole model: a, dipole structure; b, mesh; c, 3D directivity chart; d, 2D directivity chart in the plane ϕ ¼ 0.

computational region by the so-called air box with ABC on its surface. The code automatically determines air box size so that the minimal distance between the dipole and the air box is not less than λ=4, where λ is wavelength in free space.

Cylinders have diameter 0.5 mm and length 9.5 mm. Together with a lumped port of 1-mm length, the total dipole length is 20 mm. Excitation frequency was chosen to be 1.5 GHz, so the ratio l=λ ¼ 0:1. Such a dipole can be considered as elemental (Figure 6a). Its directivity is described by the formula <sup>D</sup>ðθ,ϕÞ ¼ <sup>1</sup>:5sin<sup>2</sup>θ:

Figure 6b shows mesh used in simulation. The sphere inside cylinder contains refined local mesh, surrounding the dipole. The total number of tetrahedrons is 64553. The 3D directivity chart is shown in Figure 6c, and 2D directivity chart in plain ϕ ¼ 0 in Figure 6d. Calculated dipole directivity is 1.506 with relative error of 0.4 %. This error can be caused by the finite diameter of cylinders.

8. Conclusion

Table 2. SAR, averaged over 1 g of brain tissue

Acknowledgements

Author details

References

1987;49:65–90.

Andrey D. Grigoriev

The problem of high-frequency electromagnetic field simulation is formulated. A short survey of numerical methods for solving high-frequency electromagnetic problems is presented. It was shown that one of the most efficient methods for solving inner EMP and outer EMP with moderate electrical size is the vector finite element method in frequency domain. The algorithm of this method, including mesh generation, building of the global matrix, boundary conditions approximation, SLAE solving and FFS technique was described. Some peculiarities of lumped ports and elements implementation into the RFS computer code were also depicted. Simulation results for a simple systems having analytical solution, as well as for complex problems, such as handset antenna analysis near human head show high accuracy and efficiency of the code. The code can be used for the solution of antenna problems, waveguide

900 1.08 1.06 2000 1.18 0.99

RFS SEMCAD

Computer Simulation of High‐Frequency Electromagnetic Fields

http://dx.doi.org/10.5772/67497

211

The author expresses his deep gratitude for valuable help provided by R.V. Salimov and R.I.

Saint-Petersburg State Electrotechnical University 'LETI', Saint-Petersburg, Russia

[1] Higdon R.L. Numerical absorbing boundary conditions for wave equation. Math. Comp.

problems, PCB analysis and microwave resonator problems.

Frequency, MHz SAR, W/kg

Tikhonov from the Russian R&D Centre LG Electronics Inc.

Address all correspondence to: adgrigoriev@eltech.ru

Figure 7a shows handset together with human head and hand phantoms. The handset has PIFA-type antenna. Reflection coefficient of this antenna in the presence of the phantoms was calculated by RFS code and code SEMCAD [19], based on FDTD method.

Figure 7. A handset with human head and hand phantoms: a, model; b, reflection coefficient chart.

Finite element mesh built by the RFS code in the phantom contains only 200,000 tetrahedrons. The total mesh size was about 450,000 tetrahedrons. Efficient algorithm of SAR calculation, implemented in the RFS code, provides using such comparatively small mesh size without the loss of accuracy. SEMCAD hexagonal mesh had about 1.5 billion of cells.

Figure 7b shows the module of reflection coefficient versus frequency for FRS and SEMCAD. As can be seen, both codes give similar results, especially for resonant frequencies.

We also calculated maximum SAR level in a brain tissue, averaged on 1 g of the tissue. The results are shown in Table 2. As can be seen, both codes give similar results, but SEMCAD solution time was nearly five times greater than RFS.


Table 2. SAR, averaged over 1 g of brain tissue

## 8. Conclusion

computational region by the so-called air box with ABC on its surface. The code automatically determines air box size so that the minimal distance between the dipole and the air box is not

Cylinders have diameter 0.5 mm and length 9.5 mm. Together with a lumped port of 1-mm length, the total dipole length is 20 mm. Excitation frequency was chosen to be 1.5 GHz, so the ratio l=λ ¼ 0:1. Such a dipole can be considered as elemental (Figure 6a). Its directivity is

Figure 6b shows mesh used in simulation. The sphere inside cylinder contains refined local mesh, surrounding the dipole. The total number of tetrahedrons is 64553. The 3D directivity chart is shown in Figure 6c, and 2D directivity chart in plain ϕ ¼ 0 in Figure 6d. Calculated dipole directivity is 1.506 with relative error of 0.4 %. This error can be caused by the finite

Figure 7a shows handset together with human head and hand phantoms. The handset has PIFA-type antenna. Reflection coefficient of this antenna in the presence of the phantoms was

Finite element mesh built by the RFS code in the phantom contains only 200,000 tetrahedrons. The total mesh size was about 450,000 tetrahedrons. Efficient algorithm of SAR calculation, implemented in the RFS code, provides using such comparatively small mesh size without the

Figure 7b shows the module of reflection coefficient versus frequency for FRS and SEMCAD.

We also calculated maximum SAR level in a brain tissue, averaged on 1 g of the tissue. The results are shown in Table 2. As can be seen, both codes give similar results, but SEMCAD

As can be seen, both codes give similar results, especially for resonant frequencies.

loss of accuracy. SEMCAD hexagonal mesh had about 1.5 billion of cells.

Figure 7. A handset with human head and hand phantoms: a, model; b, reflection coefficient chart.

solution time was nearly five times greater than RFS.

calculated by RFS code and code SEMCAD [19], based on FDTD method.

less than λ=4, where λ is wavelength in free space.

described by the formula <sup>D</sup>ðθ,ϕÞ ¼ <sup>1</sup>:5sin<sup>2</sup>θ:

diameter of cylinders.

210 Computer Simulation

The problem of high-frequency electromagnetic field simulation is formulated. A short survey of numerical methods for solving high-frequency electromagnetic problems is presented. It was shown that one of the most efficient methods for solving inner EMP and outer EMP with moderate electrical size is the vector finite element method in frequency domain. The algorithm of this method, including mesh generation, building of the global matrix, boundary conditions approximation, SLAE solving and FFS technique was described. Some peculiarities of lumped ports and elements implementation into the RFS computer code were also depicted. Simulation results for a simple systems having analytical solution, as well as for complex problems, such as handset antenna analysis near human head show high accuracy and efficiency of the code. The code can be used for the solution of antenna problems, waveguide problems, PCB analysis and microwave resonator problems.

## Acknowledgements

The author expresses his deep gratitude for valuable help provided by R.V. Salimov and R.I. Tikhonov from the Russian R&D Centre LG Electronics Inc.

## Author details

Andrey D. Grigoriev

Address all correspondence to: adgrigoriev@eltech.ru

Saint-Petersburg State Electrotechnical University 'LETI', Saint-Petersburg, Russia

## References

[1] Higdon R.L. Numerical absorbing boundary conditions for wave equation. Math. Comp. 1987;49:65–90.

[2] Stratton J.A. Electromagnetic Theory. 1st ed. New York and London: McGraw-Hill Book Company; 1941. 648 p.

**Chapter 10**

**Modeling and Simulation of Task Allocation with**

The task allocation problem is a key element in the solution of several applications from different engineering fields. With the explosion of the amount of information produced by the today Internet-connected solutions, scheduling techniques for the allocation of tasks relying on grids, clusters of computers, or in the cloud computing, is at the core of efficient solutions. The task allocation is an important problem within some branch of the computer sciences and operations research, where it is usually modeled as an optimization of a combinatorial problem with the inconvenience of a state explosion problem. This chapter proposes the modeling of the task allocation problem by the use of Colored Petri nets. The proposed methodology allows the construction of compact models for task scheduling problems. Moreover, a simulation process is possible within the constructed model, which allows the study of some performance aspects of the task

Keywords: colored Petri nets, task allocation problem, modeling and simulation, dis-

The services provided by Internet-based modern applications require the execution of massively parallel processes in order to produce time-effective information for different requirements. The processing of these "Big Data" is very computing, demanding in almost all of the cases. Thus, vast server farms installed over the world are ready to process the requests from different users. The digital data available in the web are huge and diverse, therefore complex,

> © 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Colored Petri Nets**

Mildreth Alcaraz-Mejia,

http://dx.doi.org/10.5772/67950

tributed computing

and in a continuous growth rate.

1. Introduction

Abstract

Raul Campos-Rodriguez and Marco Caballero-Gutierrez

Additional information is available at the end of the chapter

allocation problem before any implementation stage.


## **Modeling and Simulation of Task Allocation with Colored Petri Nets**

Mildreth Alcaraz-Mejia, Raul Campos-Rodriguez and Marco Caballero-Gutierrez

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/67950

#### Abstract

[2] Stratton J.A. Electromagnetic Theory. 1st ed. New York and London: McGraw-Hill Book

[3] Jackson J.D. Classical Electrodynamics. 3rd ed. New York: JohnWiley & Sons, Inc.; 1999. 815 p. [4] Yee K.S. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans. Antennas Propag. 1966;14(5):302–307.

[5] Taflov A., Hagness S.C. Computational Electrodynamics: The Finite-Difference Time-

[6] Weiland T. A discretization method for the solution of Maxwell's equations for sis-com-

[7] Johns P.B., Beurle R.L. Numerical solution of 2-dimensional scattering problems using a

[8] Christpoulos C. The Transmission-Line modelling Method. Oxford: Morgan and Claypool;

[9] Courant R.L. Variational methods for solution of problems of equilibrium and vibration.

[10] Jin J. The Finite Element Method in Electromagnetic. 2nd ed. New York: John Wiley &

[11] Gibson W.C. The Method of Moments in Electromagnetics. Boca Raton, FL: Taylor &

[12] Grigoriev A.D., Salimov R.V., Tikhonov R.I. Cellular handsets antenna modelling by

[13] Grigoryev A.D., Salimov R.V., Tikhonov R.I. Efficient Analysis of Full Mobile Hanset CAD Models with Automatic Correction Geometric Errors. In: 26th Annual Review of Progress in Applied Computational Electromagnetics; April 26–29; Tampere, Finland.

[14] Nedelec J.C.Mixed finite elementsinR3.Numer.Meth.1980;35:315–341.doi:10.1007/BF01396415. [15] Webb J.P. Hierarchical vector basis functions of arbitrary order for triangular and tetra-

[16] Andersen L.S., Volakis J.L. Hierarchical tangential vector finite elements for tetrahedra.

[17] Grigoryev A.D., Salimov R.V., Tikhonov R.I. Multiple-cell lumped element and port models foe the vector finite element method. Electromagnetics. 2008;25(6):18–26.

[18] Slone R.D., Lee R., Lee J.-F. Well-conditioned asymptotic wave form evaluation for finite

[19] SPEAG. SEMCAD X, intuitive, powerful, DC to light [Internet]. 2014. Retrieved from

hedral finite elements. IEEE Trans. Antennas Propag. 1999;47(8):1244–1253.

finite element method. J. Commun. Technol. Electron. 2012;53(3):261–270.

Tampere: Tampere University of Technology; 2010. pp. 416–420.

elements. IEEE Trans. Antennas Propag. 2003; 51(9):2442–2447.

IEEE Microw. Guided Wave Lett. 1998;8(3):127–129.

Domain Method. 3rd ed. New York: Artech House; 2006. 1006 p.

ponent fields. Electron. Commun. 1977;31(3):401–410.

Bull. Am. Math. Soc. 1943;5(1):1–23.

Sons, Inc.; 2002. 752 p.

Francis Group; 2008. 272 p.

http://www.speag.com/

transmission-line matrix. Proc. IIE. 1971;118(9):1203–1208.

Company; 1941. 648 p.

212 Computer Simulation

2006. 124 p.

The task allocation problem is a key element in the solution of several applications from different engineering fields. With the explosion of the amount of information produced by the today Internet-connected solutions, scheduling techniques for the allocation of tasks relying on grids, clusters of computers, or in the cloud computing, is at the core of efficient solutions. The task allocation is an important problem within some branch of the computer sciences and operations research, where it is usually modeled as an optimization of a combinatorial problem with the inconvenience of a state explosion problem. This chapter proposes the modeling of the task allocation problem by the use of Colored Petri nets. The proposed methodology allows the construction of compact models for task scheduling problems. Moreover, a simulation process is possible within the constructed model, which allows the study of some performance aspects of the task allocation problem before any implementation stage.

Keywords: colored Petri nets, task allocation problem, modeling and simulation, distributed computing

### 1. Introduction

The services provided by Internet-based modern applications require the execution of massively parallel processes in order to produce time-effective information for different requirements. The processing of these "Big Data" is very computing, demanding in almost all of the cases. Thus, vast server farms installed over the world are ready to process the requests from different users. The digital data available in the web are huge and diverse, therefore complex, and in a continuous growth rate.

© 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Nowadays, most of the queries on Internet are produced from smartphones and personal computers. However, the data traffic on web is expected to grow in an exponential rate, as the technology related to the Internet of Things (IoT) allows the connection of other devices from the houses, smart warehouses, automated driving cars, machinery from the agriculture, UAV's for goods delivering, or smart cities, to mention a few. In order to cope with the challenges related to a massive number of connected devices, new ways of computing emerged. Most of them are based on parallel and distributed computing, such as clusters, grids, and cloud computing paradigms.

graphics. The colored tokens represent the processes, tasks, and other variables of interest in the problem. A decision function is attached to some transitions of the CPN model. This function influences the behavior of the process or thread. Due to the graphical nature of the CPN, in addition to the mathematical fundamentals, this methodology provides a great framework to simulate and analyze the task allocation problem, avoiding the state space explosion, and thus, having a compact representation of a complex behavior. The flexibility of the modeling framework allows simulating different scenarios, using different decision functions, as well as varying the number of initial processes or threads. Additionally, the Petri Net model can be

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

215

The rest of this chapter includes a background of the task allocation and several techniques that have been used for addressing this problem. The modeling framework based on colored Petri nets is next presented and the modeling methodology is detailed. This methodology highlights the main aspects of the task allocation problem, such as the concepts of task and process, as well as the balancing policies and how to capture them by CPN blocks. The simulation of the obtained CPN models is then presented and discussed. Some performance charts are provided. The charts allow the study of key aspects of task allocation problem such as the effective work done by a processor, the task contention, the effects of different balancing policies, to mention a few. With this information, the scientists and engineers are able to study different aspects of the problem prior any implementation stage in a particular field of interest.

Task allocation problem has been addressed for several purposes and using different techniques. Analytic solutions are in some cases prohibitive since they require the analysis of an exponential number of possible solutions. Thus, some designs use probabilistic approaches to produce good solutions in a reasonable time. Jevtić et al. [28] present a distributed bee's algorithm for static task allocation inspired on the behavior of a colony of bees and based on the optimization of a cost function. They also show results using the Arena simulator. Delle Fave et al. [29] introduced a decentralized optimization for static task assignment based on a constrained utility function, firstly evaluated by simulations. Johnson et al. [30] state a distributed multi-agent algorithm for static multi-task allocation as an optimization problem based on mixed-integer programming. They also provide some results of Monte Carlo simulation. Zhao-Pin et al. [31] introduced a distributed and static task allocation algorithm for multiagents that learn how to determine by themselves what tasks should realize and how to form coalitions for cooperation, via their proposed profit-sharing learning mechanism. Macarthur et al. [32] introduced a dynamic and distributed algorithm for multi-agent task allocation problems, based on fast-max-sum algorithm combined with a branch-and-bound technique, in order to reduce the execution time in the allocation of task to its respective agent. Alistarh et al. [27, 33] present a decentralized task allocation, static and dynamic, based on to-do trees,

where decisions on every inner node are based on a probability function.

One of these most used techniques, due to its simplicity of implementation, is based on a rooted tree structure that the process traverses in order to acquire a task [7, 8]. The tree structure allows to the working processes for executing a non-deterministic behavior, which

easily extended to simulate a wider number of tasks.

2. Task allocation problem

At the core of these paradigms are the tasks scheduling and resources allocation problems. Adaptations of theory, methods, algorithms, tools, and others are arising in order to cope with the increasing number of cores per processor, as well as with the interconnected and heterogeneous computers, considering that, in many cases, those computing architectures are connected through the Internet. The tasks scheduling and resources allocation were well studied in the context of architectures with a relatively small number of processors. One of the fundamental concerns deals with designing algorithms for an efficiently allocation of the tasks among all the available processes and resources.

The answer to this question may lead to a hard work, since, in the worst case, it involves an exponential number of possible solutions [1]. Thus, a simulation process, for scenarios with a big number of computing elements, or cores, as in cloud and grid computing, is a valuable tool. In this context, the modeling and simulation are disciplines widely used to obtain feedback and statistics information, to improve associated applications and algorithms. Through experimentation in simulators, a fast prototyping process allows performance evaluation and parameter tuning for different workload conditions. To perform an efficient simulation process, a suitable modeling method is required. The Petri nets (PN's) and its extension, such as colored PN (CPN), have been considered as an efficient modeling technique for concurrent and distributed systems.

The simulation and implementation stages of scheduling and resource allocation problems in distributed systems have been successfully addressed [2–9]. Moreover, variants of PN's are used for modeling and simulation of other distributed concepts related to cloud services and computing [9–16], as well as in the field of the diagnosis of discrete event systems [17–19], reconfiguration [20–24], fuzzy inference [25] or identification [26].

This work proposes the modeling of the task allocation problem by the use of a methodology based on CPN's. It also shows how to simulate these models in order to obtain relevant information for the design process. To illustrate the techniques, this work considers a tree structure and supposes that the set of tasks and the set of processes are fixed and defined before the task allocation is performed as in Ref. [27]. The proposed CPN modeling method represents the behavior of a set of processors, or working threads, that traverse a tree structure from the root to a single leaf in order to acquire a task. The processes decide what route to take at each node of the tree structure, by using a decision function. Then, it routes the tree backwards updating some global variables that serves as inter-process communication mechanism.

The proposed method allows obtaining a CPN, which is a compact representation of the problem, and at the same time, it allows simulating its behavior and obtaining performance graphics. The colored tokens represent the processes, tasks, and other variables of interest in the problem. A decision function is attached to some transitions of the CPN model. This function influences the behavior of the process or thread. Due to the graphical nature of the CPN, in addition to the mathematical fundamentals, this methodology provides a great framework to simulate and analyze the task allocation problem, avoiding the state space explosion, and thus, having a compact representation of a complex behavior. The flexibility of the modeling framework allows simulating different scenarios, using different decision functions, as well as varying the number of initial processes or threads. Additionally, the Petri Net model can be easily extended to simulate a wider number of tasks.

The rest of this chapter includes a background of the task allocation and several techniques that have been used for addressing this problem. The modeling framework based on colored Petri nets is next presented and the modeling methodology is detailed. This methodology highlights the main aspects of the task allocation problem, such as the concepts of task and process, as well as the balancing policies and how to capture them by CPN blocks. The simulation of the obtained CPN models is then presented and discussed. Some performance charts are provided. The charts allow the study of key aspects of task allocation problem such as the effective work done by a processor, the task contention, the effects of different balancing policies, to mention a few. With this information, the scientists and engineers are able to study different aspects of the problem prior any implementation stage in a particular field of interest.

## 2. Task allocation problem

Nowadays, most of the queries on Internet are produced from smartphones and personal computers. However, the data traffic on web is expected to grow in an exponential rate, as the technology related to the Internet of Things (IoT) allows the connection of other devices from the houses, smart warehouses, automated driving cars, machinery from the agriculture, UAV's for goods delivering, or smart cities, to mention a few. In order to cope with the challenges related to a massive number of connected devices, new ways of computing emerged. Most of them are based on parallel and distributed computing, such as clusters,

At the core of these paradigms are the tasks scheduling and resources allocation problems. Adaptations of theory, methods, algorithms, tools, and others are arising in order to cope with the increasing number of cores per processor, as well as with the interconnected and heterogeneous computers, considering that, in many cases, those computing architectures are connected through the Internet. The tasks scheduling and resources allocation were well studied in the context of architectures with a relatively small number of processors. One of the fundamental concerns deals with designing algorithms for an efficiently allocation of the tasks among all the

The answer to this question may lead to a hard work, since, in the worst case, it involves an exponential number of possible solutions [1]. Thus, a simulation process, for scenarios with a big number of computing elements, or cores, as in cloud and grid computing, is a valuable tool. In this context, the modeling and simulation are disciplines widely used to obtain feedback and statistics information, to improve associated applications and algorithms. Through experimentation in simulators, a fast prototyping process allows performance evaluation and parameter tuning for different workload conditions. To perform an efficient simulation process, a suitable modeling method is required. The Petri nets (PN's) and its extension, such as colored PN (CPN), have been considered as an efficient modeling technique for concurrent and

The simulation and implementation stages of scheduling and resource allocation problems in distributed systems have been successfully addressed [2–9]. Moreover, variants of PN's are used for modeling and simulation of other distributed concepts related to cloud services and computing [9–16], as well as in the field of the diagnosis of discrete event systems [17–19],

This work proposes the modeling of the task allocation problem by the use of a methodology based on CPN's. It also shows how to simulate these models in order to obtain relevant information for the design process. To illustrate the techniques, this work considers a tree structure and supposes that the set of tasks and the set of processes are fixed and defined before the task allocation is performed as in Ref. [27]. The proposed CPN modeling method represents the behavior of a set of processors, or working threads, that traverse a tree structure from the root to a single leaf in order to acquire a task. The processes decide what route to take at each node of the tree structure, by using a decision function. Then, it routes the tree backwards updating some

The proposed method allows obtaining a CPN, which is a compact representation of the problem, and at the same time, it allows simulating its behavior and obtaining performance

reconfiguration [20–24], fuzzy inference [25] or identification [26].

global variables that serves as inter-process communication mechanism.

grids, and cloud computing paradigms.

214 Computer Simulation

available processes and resources.

distributed systems.

Task allocation problem has been addressed for several purposes and using different techniques. Analytic solutions are in some cases prohibitive since they require the analysis of an exponential number of possible solutions. Thus, some designs use probabilistic approaches to produce good solutions in a reasonable time. Jevtić et al. [28] present a distributed bee's algorithm for static task allocation inspired on the behavior of a colony of bees and based on the optimization of a cost function. They also show results using the Arena simulator. Delle Fave et al. [29] introduced a decentralized optimization for static task assignment based on a constrained utility function, firstly evaluated by simulations. Johnson et al. [30] state a distributed multi-agent algorithm for static multi-task allocation as an optimization problem based on mixed-integer programming. They also provide some results of Monte Carlo simulation. Zhao-Pin et al. [31] introduced a distributed and static task allocation algorithm for multiagents that learn how to determine by themselves what tasks should realize and how to form coalitions for cooperation, via their proposed profit-sharing learning mechanism. Macarthur et al. [32] introduced a dynamic and distributed algorithm for multi-agent task allocation problems, based on fast-max-sum algorithm combined with a branch-and-bound technique, in order to reduce the execution time in the allocation of task to its respective agent. Alistarh et al. [27, 33] present a decentralized task allocation, static and dynamic, based on to-do trees, where decisions on every inner node are based on a probability function.

One of these most used techniques, due to its simplicity of implementation, is based on a rooted tree structure that the process traverses in order to acquire a task [7, 8]. The tree structure allows to the working processes for executing a non-deterministic behavior, which is quite suitable for real environments. Typically, in these tree structures, every reference to a task is allocated on the leaves. Thus, every internal vertex allows to a working process taking a decision based on a function over the number of tasks in every leaf of the subtrees from each child of the current vertex. Hence, the decision function, which could be a simple deterministic algorithm, or a sophisticated statistical or stochastic procedure, is a key element for which a process traverses through the tree structure, in this kind of task allocation techniques.

left subtree and four pending tasks at the right one. Then, at every stage, the processes decide to traverse to the left or right in an asynchronous fashion. The total number of pending tasks is backward updated from the lower nodes to the uppers as the tasks are successfully executed by the threads. Notice that, the number of pending tasks at the left or right subtree is an

Since the processes are asynchronous among them, it is possible that they arrive to the leaf nodes that may occur at different rates. Every time that a process reaches a leaf, it asks for the execution of the task. If such a task is free, then the process blocks this task, executes the activity, and marks the task as executed. Then, it goes back to the root node to seek for a new task to execute, while the information of the remaining tasks is backward updated from the leaf to the root node. It is possible that when a process reaches a leaf node, the corresponding task was already executed by another process. In this case, the process goes back to the root

node, and this miss is counted to measure the effectiveness of the algorithms.

r <

x

In Eq. (1), it represented a decision rule that a working process may execute at every level of the BDT. Let x be the total number of pending tasks in the left subtree and y be the total number of pending tasks in the right subtree, while r is a random number. Thus, the probability to go to the left subtree in a flip coin is given by r. As greater the ratio, it is greater the probability to go to the left child of the current node. One of the main objectives of this work is to measure the effectiveness of a decision rule as Eq. (1) by means of a simulation process.

PN's are a modeling framework that combines a graphical visualization with a mathematical model [35]. The CPN's are one of the extensions to PN's [9, 17, 35–42], which allow the construction of compact representations of big models. In this section, basic notions of CPN

The main characteristic that differentiates CPN from other type of nets resides in the token definition. In a CPN, the tokens can stand for complex data besides of a single value, such as an integer, a real, Boolean or strings. This characteristic allows for the representation of elaborated data types, similar to those used by high-level programming languages. This ability exploits the multi-set cardinality to construct compact models that otherwise are in the power set of the

colored tokens. The formal definition of a CPN is as follows [35].

N ¼ 〈P, T, Pre, Post, B, F, Cp, Ct〉 where its elements are described as:

• T is a finite set of transitions of N , such that T∩ P ¼ ∅, with n ¼ jTj,

Definition 1. A Colored Petri Net (CPN) is a tuple

• P is a finite set of places of N , with m ¼ jPj,

<sup>x</sup> <sup>þ</sup> <sup>y</sup> <sup>ð</sup>1<sup>Þ</sup>

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

217

approximation since this number may be not updated yet.

3. Colored Petri nets

are presented.

Figure 1 shows a binary tree that represents a scheduling problem where eight tasks have to be attended. The height of the tree is H ¼ log2ð8Þ ¼ 3, where the level of v2, <sup>2</sup> is 2, the vertex v1,<sup>0</sup> is said to be parent, or ancestor, of v2,<sup>0</sup> and v2, 1, and the latter two are said to be children, or descendant, of v1,0. The vertex v0,<sup>0</sup> is the root, and every vertex is internal except for all of the level 3, i.e., v3, <sup>x</sup> with x∈½0, 7�, which are known as leafs of the tree. The tree formed by the vertexes v2, <sup>1</sup>, v3, <sup>2</sup>, v3, <sup>3</sup> and the corresponding edges that connects them, is a subtree with v2,<sup>1</sup> as a root.

From the graph theory point of view, a DT is an acyclic graph in which a unique path connects every node to others, represented in a top-to-bottom fashion, where the root, i.e., initial node, is plotted at the top of the tree. The DT is a binary decision tree (BDT) if at every node is considered at most two possible decisions, as shown in Figure 1. For more information on trees see Ref. [34].

#### 2.1. Task scheduling with DT

By using a tree structure like the one depicted in Figure 1, the task allocation algorithm is relatively simple. Each working process shares the tree structure and iterates in an infinite loop, trying to acquire a new task to execute, which is located at the leafs of the tree. Each working process starts at the root node of the tree and decides whether to go to the left or the right child node. The process knows an approximation of the total pending tasks at the left and right subtrees at every step. For example, in Figure 1, each working thread knows that at the beginning of the scheduling process, at the root node, that there exist four pending tasks at the

Figure 1. A direct rooted full binary tree.

left subtree and four pending tasks at the right one. Then, at every stage, the processes decide to traverse to the left or right in an asynchronous fashion. The total number of pending tasks is backward updated from the lower nodes to the uppers as the tasks are successfully executed by the threads. Notice that, the number of pending tasks at the left or right subtree is an approximation since this number may be not updated yet.

Since the processes are asynchronous among them, it is possible that they arrive to the leaf nodes that may occur at different rates. Every time that a process reaches a leaf, it asks for the execution of the task. If such a task is free, then the process blocks this task, executes the activity, and marks the task as executed. Then, it goes back to the root node to seek for a new task to execute, while the information of the remaining tasks is backward updated from the leaf to the root node. It is possible that when a process reaches a leaf node, the corresponding task was already executed by another process. In this case, the process goes back to the root node, and this miss is counted to measure the effectiveness of the algorithms.

$$r < \frac{x}{x+y} \tag{1}$$

In Eq. (1), it represented a decision rule that a working process may execute at every level of the BDT. Let x be the total number of pending tasks in the left subtree and y be the total number of pending tasks in the right subtree, while r is a random number. Thus, the probability to go to the left subtree in a flip coin is given by r. As greater the ratio, it is greater the probability to go to the left child of the current node. One of the main objectives of this work is to measure the effectiveness of a decision rule as Eq. (1) by means of a simulation process.

## 3. Colored Petri nets

is quite suitable for real environments. Typically, in these tree structures, every reference to a task is allocated on the leaves. Thus, every internal vertex allows to a working process taking a decision based on a function over the number of tasks in every leaf of the subtrees from each child of the current vertex. Hence, the decision function, which could be a simple deterministic algorithm, or a sophisticated statistical or stochastic procedure, is a key element for which a

Figure 1 shows a binary tree that represents a scheduling problem where eight tasks have to be attended. The height of the tree is H ¼ log2ð8Þ ¼ 3, where the level of v2, <sup>2</sup> is 2, the vertex v1,<sup>0</sup> is said to be parent, or ancestor, of v2,<sup>0</sup> and v2, 1, and the latter two are said to be children, or descendant, of v1,0. The vertex v0,<sup>0</sup> is the root, and every vertex is internal except for all of the level 3, i.e., v3, <sup>x</sup> with x∈½0, 7�, which are known as leafs of the tree. The tree formed by the vertexes v2, <sup>1</sup>, v3, <sup>2</sup>, v3, <sup>3</sup> and the corresponding edges that connects them, is a subtree with v2,<sup>1</sup> as

From the graph theory point of view, a DT is an acyclic graph in which a unique path connects every node to others, represented in a top-to-bottom fashion, where the root, i.e., initial node, is plotted at the top of the tree. The DT is a binary decision tree (BDT) if at every node is considered at most two possible decisions, as shown in Figure 1. For more information on trees

By using a tree structure like the one depicted in Figure 1, the task allocation algorithm is relatively simple. Each working process shares the tree structure and iterates in an infinite loop, trying to acquire a new task to execute, which is located at the leafs of the tree. Each working process starts at the root node of the tree and decides whether to go to the left or the right child node. The process knows an approximation of the total pending tasks at the left and right subtrees at every step. For example, in Figure 1, each working thread knows that at the beginning of the scheduling process, at the root node, that there exist four pending tasks at the

process traverses through the tree structure, in this kind of task allocation techniques.

a root.

216 Computer Simulation

see Ref. [34].

2.1. Task scheduling with DT

Figure 1. A direct rooted full binary tree.

PN's are a modeling framework that combines a graphical visualization with a mathematical model [35]. The CPN's are one of the extensions to PN's [9, 17, 35–42], which allow the construction of compact representations of big models. In this section, basic notions of CPN are presented.

The main characteristic that differentiates CPN from other type of nets resides in the token definition. In a CPN, the tokens can stand for complex data besides of a single value, such as an integer, a real, Boolean or strings. This characteristic allows for the representation of elaborated data types, similar to those used by high-level programming languages. This ability exploits the multi-set cardinality to construct compact models that otherwise are in the power set of the colored tokens. The formal definition of a CPN is as follows [35].

Definition 1. A Colored Petri Net (CPN) is a tuple

N ¼ 〈P, T, Pre, Post, B, F, Cp, Ct〉 where its elements are described as:


• Cpðp1Þ ¼ Cpðp2Þ ¼ o, Cpðp3Þ ¼ r,

• Pre and Post are as depicted, and

0

and <sup>M</sup><sup>1</sup>! <sup>t</sup>2,ðo2<sup>Þ</sup>

for the color condition <sup>ð</sup>o1, r2Þ, then <sup>M</sup><sup>1</sup> ! <sup>t</sup>1,ðo1,r2<sup>Þ</sup>

function intended for the characters in this CPN model.

4. Modeling and simulation of task scheduling

4.1. Structure of the task allocation as a CPN

, where o<sup>1</sup> ¼ 1, o<sup>2</sup> ¼ 2, r<sup>1</sup> ¼ a, r<sup>2</sup> ¼ b.

∅, {r1,r2}� : {ðo1,r1Þðo2,r2Þ}, thus, t<sup>1</sup> can be fired for any color binding fðo1, r2Þðo2, r1Þg due to the

Roughly speaking, the CPN binds the even number in p<sup>1</sup> to the character "a" in p<sup>3</sup> and the odd number to character "b." This is defined by the guard function f <sup>1</sup> attached to t1. The function s<sup>1</sup> discards the r element of the token bound to the firing of t1, while the firing of t<sup>2</sup> adds two to the o element of the token and recovers the r element that is put back to the place p3. Thus, the tones in p<sup>1</sup> changes to o<sup>1</sup> ¼ 1, 3, 5, … and o<sup>2</sup> ¼ 2, 4, 6, …, while the characters "a" and "b" in p<sup>3</sup> remain the same all the time. The "R" in the place p<sup>3</sup> stands for "resources," which is the

Notice that, the bindings ðo1, r1Þ and ðo2, r2Þ are discarded by the function f <sup>1</sup>. This is one of the advantages of the CPN modeling tool, since otherwise, the combinatorial explosion of possible bindings turns untreatable under some circumstances that typically arises in real applications. For more information about the CPN modeling formalism and related analysis and tools, see

The construction of a generalized model for a BDT based on CPN considers d different features and nT classes. The procedure includes the identification of the transitions and places of the model, the arcs and its labeling, as well as the multi-set of colored tokens. The remaining of this

Consider the following steps for the construction of a CPN model that captures the structure of a BDT. It is supposed that the BDT is a full binary tree with a complete set of tasks represented by a structure called Task Array. Thus, the size of the Task Array is 2<sup>H</sup>, where the height of the

section is devoted to detail the methodology and illustrate it by short examples.

, where <sup>o</sup><sup>1</sup> <sup>¼</sup> 1, <sup>o</sup><sup>2</sup> <sup>¼</sup> 2, <sup>r</sup><sup>2</sup> <sup>¼</sup> <sup>b</sup>, by Post½■, t1�ð<sup>f</sup> <sup>1</sup>Þ ¼ {ðs1Þ} as detailed in the

, where o<sup>1</sup> ¼ 1, o<sup>2</sup> ¼ 4, r<sup>1</sup> ¼ a, r<sup>2</sup> ¼ b. However, if t<sup>1</sup> is fired from M<sup>1</sup>

, in such a way, that if <sup>t</sup><sup>2</sup> is fired, then <sup>M</sup><sup>1</sup> !<sup>t</sup>2,ðo2<sup>Þ</sup>

M3, with M<sup>3</sup> ¼ ½∅, fo1, o2g,∅�

Modeling and Simulation of Task Allocation with Colored Petri Nets

≥ ½fo1, o2g,

219

M1, with

M2, with

, where o<sup>1</sup> ¼ 1,

0

http://dx.doi.org/10.5772/67950

Notice that, at the initial condition, <sup>M</sup><sup>0</sup> <sup>≥</sup> Pre½■, t1�ð<sup>f</sup> <sup>1</sup>Þ. That is ½fo1, o2g, <sup>∅</sup>, <sup>f</sup>r1, r2g�<sup>0</sup>

condition in <sup>f</sup>1. Suppose that, <sup>t</sup><sup>1</sup> fires for color condition <sup>f</sup><sup>1</sup> with <sup>ð</sup>o2, r1Þ, then <sup>M</sup><sup>0</sup> ! <sup>t</sup>1,ðo2,r1<sup>Þ</sup>

• Ctðt1Þ ¼ f <sup>1</sup>, Ctðt2Þ ¼ f <sup>2</sup>,

• M<sup>0</sup> ¼ ½fo1, o2g, ∅, {r1,r2}�

M<sup>1</sup> ¼ ½fo1g, fo2g, fr2g�<sup>0</sup>

figure. Now, <sup>M</sup><sup>1</sup> ! <sup>t</sup>1,ðo1,r2<sup>Þ</sup>

<sup>M</sup><sup>2</sup> ¼ ½fo1, <sup>o</sup>2g,∅,fr1,r2g�<sup>0</sup>

o<sup>2</sup> ¼ 2.

Ref. [36].

tree is H.

The incidence matrix is C ¼ Post � Pre. The mapping given by Pre½p, t� : CtðtÞ ! BagðCpðpÞÞ, defines for each conditional color mapping f of t ði:e: ∀f <sup>i</sup> ∈ CtðtÞÞ, the token bag to be removed from p, in the occurrence of t, under color condition f. In the same way, Post½p, t�ðfÞ specifies the token bag to be added to p, in the occurrence of t, under color condition f.


Definition 4. The marking evolution w.r.t a color condition <sup>f</sup> is given by <sup>M</sup><sup>0</sup> <sup>¼</sup> <sup>M</sup> <sup>þ</sup> Post½■, t�ðf<sup>Þ</sup>

�Pre½■, t�ðfÞ ¼ <sup>M</sup> <sup>þ</sup> <sup>C</sup>½■, t�ðfÞ, denoted by <sup>M</sup>! t, f M<sup>0</sup> . For a general color condition, it is denoted as M! t M<sup>0</sup> where f is implicit in a context.

The net in Figure 2 represents a very simple CPN System ðN , M0Þ where:


Figure 2. A CPN example.


• B is finite set of color classes,

is known as a CPN System.

t

denoted as M!

• P ¼ fp1, p2, p3g,

• T ¼ ft1, t2g,

• F ¼ {f <sup>1</sup>, f <sup>2</sup>},

Figure 2. A CPN example.

• Cp : P ! B is the color domain mapping,

• Ct : T ! F is a conditional color mapping, and

in a marking <sup>M</sup> iff <sup>M</sup> <sup>≥</sup> Pre½■, t�ðfÞ, denoted as <sup>M</sup>!

�Pre½■, t�ðfÞ ¼ <sup>M</sup> <sup>þ</sup> <sup>C</sup>½■, t�ðfÞ, denoted by <sup>M</sup>!

• B ¼ o∪ r, where o ¼ fo1, o2g and r ¼ fr1, r2g are variables,

• Pre, Post∈SjP<sup>j</sup> · <sup>j</sup>T<sup>j</sup> are matrices representing the input and output incidence matrices of N , such that Pre½p, t� : CtðtÞ ! BagðCpðpÞÞ and Post½p, t� : CtðtÞ ! BagðCpðpÞÞ are mappings for every pair <sup>ð</sup>p, tÞ<sup>∈</sup> <sup>P</sup> · <sup>T</sup>; where BagðSÞ ¼ <sup>∪</sup> <sup>b</sup>ðsi<sup>Þ</sup> for all si <sup>∈</sup>S, such that <sup>b</sup>ðsiÞ ¼ <sup>X</sup>si <sup>∈</sup>S.

The incidence matrix is C ¼ Post � Pre. The mapping given by Pre½p, t� : CtðtÞ ! BagðCpðpÞÞ, defines for each conditional color mapping f of t ði:e: ∀f <sup>i</sup> ∈ CtðtÞÞ, the token bag to be removed from p, in the occurrence of t, under color condition f. In the same way, Post½p, t�ðfÞ specifies the

Definition 2. A marking M of a CPN N is a vector of size m ¼ jPj, such that MðpÞ∈ BagðCpðpÞÞ for every p ∈P. The vector M<sup>0</sup> denotes the initial marking of the net N and the pair ðN , M0Þ

Definition 3. A transition in a CPN System ðN , M0Þ is said to be enabled for a color condition f

Definition 4. The marking evolution w.r.t a color condition <sup>f</sup> is given by <sup>M</sup><sup>0</sup> <sup>¼</sup> <sup>M</sup> <sup>þ</sup> Post½■, t�ðf<sup>Þ</sup>

t, f .

t, f M<sup>0</sup>

. For a general color condition, it is

token bag to be added to p, in the occurrence of t, under color condition f.

M<sup>0</sup> where f is implicit in a context.

The net in Figure 2 represents a very simple CPN System ðN , M0Þ where:

• F is a set of conditions,

218 Computer Simulation


Notice that, at the initial condition, <sup>M</sup><sup>0</sup> <sup>≥</sup> Pre½■, t1�ð<sup>f</sup> <sup>1</sup>Þ. That is ½fo1, o2g, <sup>∅</sup>, <sup>f</sup>r1, r2g�<sup>0</sup> ≥ ½fo1, o2g, ∅, {r1,r2}� : {ðo1,r1Þðo2,r2Þ}, thus, t<sup>1</sup> can be fired for any color binding fðo1, r2Þðo2, r1Þg due to the condition in <sup>f</sup>1. Suppose that, <sup>t</sup><sup>1</sup> fires for color condition <sup>f</sup><sup>1</sup> with <sup>ð</sup>o2, r1Þ, then <sup>M</sup><sup>0</sup> ! <sup>t</sup>1,ðo2,r1<sup>Þ</sup> M1, with M<sup>1</sup> ¼ ½fo1g, fo2g, fr2g�<sup>0</sup> , where <sup>o</sup><sup>1</sup> <sup>¼</sup> 1, <sup>o</sup><sup>2</sup> <sup>¼</sup> 2, <sup>r</sup><sup>2</sup> <sup>¼</sup> <sup>b</sup>, by Post½■, t1�ð<sup>f</sup> <sup>1</sup>Þ ¼ {ðs1Þ} as detailed in the figure. Now, <sup>M</sup><sup>1</sup> ! <sup>t</sup>1,ðo1,r2<sup>Þ</sup> and <sup>M</sup><sup>1</sup>! <sup>t</sup>2,ðo2<sup>Þ</sup> , in such a way, that if <sup>t</sup><sup>2</sup> is fired, then <sup>M</sup><sup>1</sup> !<sup>t</sup>2,ðo2<sup>Þ</sup> M2, with <sup>M</sup><sup>2</sup> ¼ ½fo1, <sup>o</sup>2g,∅,fr1,r2g�<sup>0</sup> , where o<sup>1</sup> ¼ 1, o<sup>2</sup> ¼ 4, r<sup>1</sup> ¼ a, r<sup>2</sup> ¼ b. However, if t<sup>1</sup> is fired from M<sup>1</sup> for the color condition <sup>ð</sup>o1, r2Þ, then <sup>M</sup><sup>1</sup> ! <sup>t</sup>1,ðo1,r2<sup>Þ</sup> M3, with M<sup>3</sup> ¼ ½∅, fo1, o2g,∅� 0 , where o<sup>1</sup> ¼ 1, o<sup>2</sup> ¼ 2.

Roughly speaking, the CPN binds the even number in p<sup>1</sup> to the character "a" in p<sup>3</sup> and the odd number to character "b." This is defined by the guard function f <sup>1</sup> attached to t1. The function s<sup>1</sup> discards the r element of the token bound to the firing of t1, while the firing of t<sup>2</sup> adds two to the o element of the token and recovers the r element that is put back to the place p3. Thus, the tones in p<sup>1</sup> changes to o<sup>1</sup> ¼ 1, 3, 5, … and o<sup>2</sup> ¼ 2, 4, 6, …, while the characters "a" and "b" in p<sup>3</sup> remain the same all the time. The "R" in the place p<sup>3</sup> stands for "resources," which is the function intended for the characters in this CPN model.

Notice that, the bindings ðo1, r1Þ and ðo2, r2Þ are discarded by the function f <sup>1</sup>. This is one of the advantages of the CPN modeling tool, since otherwise, the combinatorial explosion of possible bindings turns untreatable under some circumstances that typically arises in real applications. For more information about the CPN modeling formalism and related analysis and tools, see Ref. [36].

## 4. Modeling and simulation of task scheduling

The construction of a generalized model for a BDT based on CPN considers d different features and nT classes. The procedure includes the identification of the transitions and places of the model, the arcs and its labeling, as well as the multi-set of colored tokens. The remaining of this section is devoted to detail the methodology and illustrate it by short examples.

#### 4.1. Structure of the task allocation as a CPN

Consider the following steps for the construction of a CPN model that captures the structure of a BDT. It is supposed that the BDT is a full binary tree with a complete set of tasks represented by a structure called Task Array. Thus, the size of the Task Array is 2<sup>H</sup>, where the height of the tree is H.

## Step 1. Labeling nodes in the BDT

Firstly, suppose that the DT is a full binary tree. Then, assign "0" to every left edge and "1" to every right edge. This labeling represents the result of the flip coin in the task allocation procedure of a BDT discussed in the previous section.

• B ¼ o ∪ r ∪ a, for o ¼ ∪ oi, r ¼ ∪ rj, and a ¼ ∪ al where oi ¼ 〈depthi, posi, vali, successi〉, rj ¼

• F ¼ ff <sup>1</sup>, f <sup>2</sup>, f <sup>3</sup>, f <sup>4</sup>g, such that f <sup>1</sup> : ðoi, rj, rkÞjdepthi þ 1 ¼ depthj ¼ depthk ∧ posi�2 ¼ posj ∧ posi�2

• CtðtiÞ ¼ f <sup>1</sup> for i∈½1, H�, CtðtiÞ ¼ f <sup>2</sup> for i∈ ½H þ 2, jPOj � 1�, CtðtiÞ ¼ f <sup>3</sup> for i ¼ H þ 1 and i ¼

Þ ¼ r for all i∈½jPOj þ 1, jPj� and Cpðpi

, ti�ðf <sup>3</sup>Þ ¼ si1jf <sup>4</sup>, for i ¼ jPOj þ 1, Pre½pi

, ti�ð<sup>f</sup> <sup>1</sup>Þ ¼ <sup>2</sup>ri1j<sup>f</sup> <sup>1</sup>, for <sup>i</sup><sup>∈</sup> <sup>½</sup>1, H�, Pre½puþ<sup>i</sup>

, ti�ð<sup>f</sup> <sup>1</sup>Þ ¼ ro1, for <sup>i</sup><sup>∈</sup> <sup>½</sup>1, H�, Post½puþ<sup>i</sup>

þ1 ¼ posk, f <sup>2</sup> : ðoi, rjÞjdepthi ¼ depthj ∧ posi ¼ posj, f <sup>3</sup> : ðoi, ajÞjposi ¼ posj, f <sup>4</sup> : ðojÞ.

• si1: o ! o, such that si1ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi, posi, vali, successi〉.

• ri1: r · r ! r, such that ri1ð〈depthj, posj, dataj〉,〈depthk, posk, datak〉Þ ¼ 〈depthj, posj, dataj〉

• Post½piþ<sup>1</sup>, ti�ð<sup>f</sup> <sup>1</sup>Þ ¼ so1j<sup>f</sup> <sup>1</sup>, for <sup>i</sup> <sup>∈</sup>½1, H�, Post½piþ<sup>1</sup>, ti�ð<sup>f</sup> <sup>2</sup>Þ ¼ so2j<sup>f</sup> <sup>2</sup>, for <sup>i</sup> <sup>∈</sup>½<sup>H</sup> <sup>þ</sup> <sup>2</sup>, <sup>j</sup>POj � <sup>1</sup>�, Post

• so<sup>1</sup> : o ! o, such that so1ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi þ 1, posi�2 þ hi, vali þ ðhiÞ�

, successi〉, withhi ¼ 0 if αholds, otherwise 1. α is a decision function.

• so<sup>3</sup> : o ! o, such that so3ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi, posi, vali, successi ¼ 1if Exist ðposi, alÞ ¼ 1 otherwise 0〉, where Existðposi, alÞ returns 1 if for a color class al ¼ posl, existl

• qo1: a ! a, such that qo1ð〈posl, existl〉Þ¼ð〈posl, existl ¼ 0 if Existðposl, alÞ ¼ 1, otherwise1〉Þ: • ro1: r · r ! r, such that ro1ð〈depthj, posj, dataj〉,〈depthk, posk, datak〉Þ ¼ 〈depthj, posj, dataj〉

• ro2: r ! r, such that ro2ð〈depthj, posj, dataj〉Þ ¼ 〈depthj, posj, dataj � 1 if successi ¼ 1,

posi, vali, successi, dataj, existl variables and depthj, posj, posl constants.

〉 for <sup>i</sup><sup>∈</sup> <sup>½</sup>1, d�, <sup>j</sup><sup>∈</sup> <sup>½</sup>1, <sup>2</sup><sup>H</sup>þ<sup>1</sup> � <sup>2</sup>�, and <sup>l</sup><sup>∈</sup> <sup>½</sup>1, nT�, with depthi,

Modeling and Simulation of Task Allocation with Colored Petri Nets

, ti�ðf <sup>2</sup>Þ ¼ si1jf <sup>2</sup>, for i∈½H þ 2, jPOj � 1�, Pre½pi

Þ ¼ a for i ¼ 2ðHþ

http://dx.doi.org/10.5772/67950

, tj�ðf <sup>3</sup>Þ ¼ qi1jf <sup>3</sup>, for

, tu�<sup>i</sup>�ðf <sup>2</sup>Þ ¼ 2ri<sup>2</sup>

, tj�ðf <sup>3</sup>Þ ¼ qo1jf <sup>3</sup>,

, tu�<sup>i</sup>�ðf <sup>2</sup>Þ ¼ ro2,

<sup>2</sup> , vali, successi〉.

, ti�ðf <sup>3</sup>Þ

221

〈depthj, posj, dataj〉, al ¼ 〈posl, readyl

Þ ¼ o for all i∈ ½1, jPOj�, Cpðpi

, ti�ðf <sup>1</sup>Þ ¼ si1jf <sup>1</sup>, for i ∈½1, H�, Pre½pi

• qi1: a ! a, such that i1ð〈posl, existlÞ〉 ¼ ð〈posl, existl〉Þ.

• ri<sup>2</sup> : r ! r, such that ri2ð〈depthj, posj, dataj〉Þ ¼ 〈depthj, posj, dataj〉.

<sup>½</sup>piþ<sup>1</sup>, ti�ð<sup>f</sup> <sup>3</sup>Þ ¼ so3j<sup>f</sup> <sup>3</sup>, for <sup>i</sup> <sup>¼</sup> <sup>H</sup> <sup>þ</sup> 1, Post½p1, ti�ð<sup>f</sup> <sup>3</sup>Þ ¼ so4j<sup>f</sup> <sup>4</sup>, for <sup>i</sup> ¼ jPOj, Post½pi

• so<sup>2</sup> : <sup>o</sup> ! <sup>o</sup>, such that so2ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi � <sup>1</sup>, posi

• so<sup>4</sup> : o ! o, such that so4ð〈depthi, posi, vali, successi〉Þ ¼ 〈0, 0, 0, 0〉,

¼ si1jf <sup>3</sup>, for i ¼ H þ 1, Pre½pi

jf <sup>2</sup>, for i ∈½1, H�, where:

〈depthk, posk, datak〉.

for i∈ ½1, H�, where:

〈depthk, posk, datak〉.

otherwisedataj〉.

2<sup>H</sup>�depthiþ<sup>1</sup>

<sup>i</sup> ¼ jPOjþjPRj þ 1, <sup>j</sup> <sup>¼</sup> <sup>H</sup> <sup>þ</sup> 1, Pre½puþ<sup>i</sup>

for <sup>i</sup> ¼ jPOjþjPRj þ 1, <sup>j</sup> <sup>¼</sup> <sup>H</sup> <sup>þ</sup> <sup>1</sup>Post½puþ<sup>i</sup>

with posl ¼ posi, existl ¼ 1, otherwise 0.

• Cpðpi

• Pre½pi

1Þ þ H þ 1.

jPOj þ 1:

After that, some information is used to label every node with a pair ðdepth, positionÞ, where depth is the level of the node, and position is the corresponding position from left to right. For example, consider the root node ð0, 0Þ in the section of the BDT depicted in Figure 3. It is a fraction of the tree of height H ¼ 3 of Figure 1. Then, to label the children of the root node proceed as follows:

• Let be depth ¼ depth � of � the � parent þ 1 and position ¼ 2 � position � of � the � parentþ h, with h ¼ 0 if this node is the left child, h ¼ 1 otherwise.

Therefore, the label of the left child is ðdepth ¼ 0 þ 1, position ¼ 0 � 2 þ 0Þ¼ð1, 0Þ, while the label of the right child is ðdepth ¼ 0 þ 1, position ¼ 0 � 2 þ 1Þ¼ð1, 1Þ.

#### Step 2. Building the CPN System

Let ðN , M0Þ be a CPN System for a dynamic binary decision tree for nT classes and d features with N ¼ 〈P, T, Pre, Post, B, F, Cp, Ct〉 constructed as follows:

• P ¼ PO∪ PR∪ PA, where PO is the set of places with assigned tokens of type O, PR is the set of places with assigned tokens of type R, PA is the place with assigned tokens of type A. jPOj ¼ 2ðH þ 1Þ, jPRj ¼ H and jPAj ¼ 1, where H ¼ log2nT. Places of type O represent every level from the root to the leaves, forward, and backwards, i.e., a token in one of those place contains the information of pos, depth, val. Places of type R represent the state of every subtree from every node at that level, i.e., the value of the data function for every node of the same level given by pos and depth. Place of type A represent the Task Array state, i.e., the information of the availability of the task.

• T is the set of transitions with jTj¼jPOj.

Figure 3. Relation BDT and CPN.

	- si1: o ! o, such that si1ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi, posi, vali, successi〉.
	- qi1: a ! a, such that i1ð〈posl, existlÞ〉 ¼ ð〈posl, existl〉Þ.

Step 1. Labeling nodes in the BDT

220 Computer Simulation

Step 2. Building the CPN System

procedure of a BDT discussed in the previous section.

with h ¼ 0 if this node is the left child, h ¼ 1 otherwise.

with N ¼ 〈P, T, Pre, Post, B, F, Cp, Ct〉 constructed as follows:

information of the availability of the task.

• T is the set of transitions with jTj¼jPOj.

Figure 3. Relation BDT and CPN.

label of the right child is ðdepth ¼ 0 þ 1, position ¼ 0 � 2 þ 1Þ¼ð1, 1Þ.

Firstly, suppose that the DT is a full binary tree. Then, assign "0" to every left edge and "1" to every right edge. This labeling represents the result of the flip coin in the task allocation

After that, some information is used to label every node with a pair ðdepth, positionÞ, where depth is the level of the node, and position is the corresponding position from left to right. For example, consider the root node ð0, 0Þ in the section of the BDT depicted in Figure 3. It is a fraction of the tree of height H ¼ 3 of Figure 1. Then, to label the children of the root node proceed as follows: • Let be depth ¼ depth � of � the � parent þ 1 and position ¼ 2 � position � of � the � parentþ h,

Therefore, the label of the left child is ðdepth ¼ 0 þ 1, position ¼ 0 � 2 þ 0Þ¼ð1, 0Þ, while the

Let ðN , M0Þ be a CPN System for a dynamic binary decision tree for nT classes and d features

• P ¼ PO∪ PR∪ PA, where PO is the set of places with assigned tokens of type O, PR is the set of places with assigned tokens of type R, PA is the place with assigned tokens of type A. jPOj ¼ 2ðH þ 1Þ, jPRj ¼ H and jPAj ¼ 1, where H ¼ log2nT. Places of type O represent every level from the root to the leaves, forward, and backwards, i.e., a token in one of those place contains the information of pos, depth, val. Places of type R represent the state of every subtree from every node at that level, i.e., the value of the data function for every node of the same level given by pos and depth. Place of type A represent the Task Array state, i.e., the

	- so<sup>1</sup> : o ! o, such that so1ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi þ 1, posi�2 þ hi, vali þ ðhiÞ� 2<sup>H</sup>�depthiþ<sup>1</sup> , successi〉, withhi ¼ 0 if αholds, otherwise 1. α is a decision function.
	- so<sup>2</sup> : <sup>o</sup> ! <sup>o</sup>, such that so2ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi � <sup>1</sup>, posi <sup>2</sup> , vali, successi〉.
	- so<sup>3</sup> : o ! o, such that so3ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi, posi, vali, successi ¼ 1if Exist ðposi, alÞ ¼ 1 otherwise 0〉, where Existðposi, alÞ returns 1 if for a color class al ¼ posl, existl with posl ¼ posi, existl ¼ 1, otherwise 0.
	- so<sup>4</sup> : o ! o, such that so4ð〈depthi, posi, vali, successi〉Þ ¼ 〈0, 0, 0, 0〉,
	- qo1: a ! a, such that qo1ð〈posl, existl〉Þ¼ð〈posl, existl ¼ 0 if Existðposl, alÞ ¼ 1, otherwise1〉Þ:
	- ro1: r · r ! r, such that ro1ð〈depthj, posj, dataj〉,〈depthk, posk, datak〉Þ ¼ 〈depthj, posj, dataj〉 〈depthk, posk, datak〉.
	- ro2: r ! r, such that ro2ð〈depthj, posj, dataj〉Þ ¼ 〈depthj, posj, dataj � 1 if successi ¼ 1, otherwisedataj〉.
	- M0½1�¼fo1, …, odg, where oi ¼ depthi ¼ 0, posi ¼ 0, vali ¼ 0, successi ¼ 0 ¼ 0, 0, 0, 0 for i∈ ½1, d�.

here introduced, it is assumed that there exists a data structure called Task Array, which may have a reference to a task, if it is initially inserted, and the task has not been taken yet. This Task Array is mapped to the location of the leaves of a tree of height H from left to right, i.e., the array has a capacity for referencing a maximum of The quaternion semi � widel tasks. Figure 1 shows a tree of height 3, with 8 leaves that are associated with the location in the Task Array, where references to task are T0, T1…, T7. Accordingly, the following functions are defined.

Definition 8. The data function is defined from the set of vertexes V of a tree T to the set of natural numbers N as data : V ! N, where dataðvÞ ¼ # of leaves in the subtree of root v with

Definition 9. The parent is a function defined from the set of vertexes of a tree T to itself as

Definition 10. The val is a recursive function defined from the set of vertexes V of a tree T of

In Figure 2, for example, valðvÞ is shown in the left of every vertex v as valðvÞ=, e.g.,

Observe that function val generates the value of a node based on the value of his parent. Notice that, there are 8 leaves in level 3, and then, every node value is exactly one entry in the Task Array. For example, the value 4= at v3, <sup>5</sup> points to the entry T4. Now, consider the nodes in level 2, which are 4, exactly the half of those at level 3. Every vertex at this level, points to the first entry in a range of 2, i.e., v2, <sup>0</sup> has a value of 0, v2, <sup>1</sup> and has a value of 2. Thus, every vertex at level 1, points to the first entry of the two partitions of the Task Array, and therefore, v0, <sup>0</sup> points

For illustration purposes, suppose that there is one node with two threads and eight tasks to be executed. Accordingly, d ¼ 2 and nT ¼ 8. The CPN System ðN , M0Þ for modeling this problem

Notice that accordingly to the proposed method, a full binary tree is obtained when nT <sup>¼</sup> <sup>2</sup><sup>n</sup> for some n, since, in this case, H ¼ log2nT is exactly equal to log2nT. In the particular example, <sup>H</sup> <sup>¼</sup> 3, since nT <sup>¼</sup> 23 <sup>¼</sup> 8. Then, <sup>j</sup>POj ¼ <sup>2</sup>ð<sup>H</sup> <sup>þ</sup> <sup>1</sup>Þ ¼ 8, with PO ¼ fp1, p2, …, p8g; <sup>j</sup>PRj ¼ <sup>H</sup> <sup>¼</sup> 3, with PR ¼ fp9, p10, p11g; jPAj ¼ 1 with PA ¼ {p12}; jTj¼jPOj ¼ 8, with T ¼ {t1, t2, …, t8}; o ¼ {o1, <sup>o</sup>2} since <sup>d</sup> <sup>¼</sup> 2; <sup>r</sup> ¼ fr1, r2, …, r14g, since 2<sup>H</sup>þ<sup>1</sup> � <sup>2</sup> <sup>¼</sup> 23 � <sup>2</sup> <sup>¼</sup> 14, and <sup>a</sup> <sup>¼</sup> {a1, <sup>a</sup>2, …, <sup>a</sup>8}.

eposðparentðv2,3Þ, depthðv2,3Þ, <sup>3</sup>

eposðparentðvi,jÞ, depthðvi,jÞ,j

�23�<sup>1</sup> <sup>þ</sup> <sup>w</sup>

 �2<sup>H</sup>�<sup>i</sup>

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

223

�23�<sup>2</sup> <sup>¼</sup> valðv1, <sup>1</sup>Þ þ <sup>w</sup>

eposðv1,1Þ, <sup>2</sup>,<sup>3</sup>

, for every pair ði, jÞ,

eposðv1, <sup>1</sup>Þ,2,<sup>3</sup>

�2<sup>1</sup> <sup>¼</sup> valðv0,0Þþ¼ <sup>0</sup>

 �21

a task ready to be taken.

• (Base): valðv0, <sup>0</sup>Þ ¼ 0.

parentðv1,1Þ

is shown in Figure 4.

valðv2,3Þ ¼ val

<sup>¼</sup> val

• (Induction): valðvi,jÞ ¼ val

where 1 ≤ i ≤ H, 0 ≤ j < 2<sup>i</sup>

parentðv2,3Þ

þwðe0,1,1Þ�4 þ wðe1,2, <sup>3</sup>Þ�2 ¼ 0 þ 1�4 þ 1�2 ¼ 6.

to the first entry of the whole of the Task Array.

The following functions complement the CPN model:

 þ w 

parent : V ! V, such that parentðvÞ ¼ y, if v is child of y.

.

 þ w 

height H, to the set of natural numbers N as val : V ! N, such that:

parentðvi,jÞ

eposðparentðv1,1Þ, depthðv1,1Þ, <sup>1</sup>

 þ w 


Figure 4 shows the relation of the tree structure with the proposed CPN as stated by the previous definition. Observe that places of type O represent every level from the root to the leaves, forward, and backwards, i.e., a token in one of those place contains the information of pos, depth, and val. Places of type R represent the state of every subtree from every node at that level, i.e., the value of the data function for every node of the same level given by pos and depth. Place of type A represents the Task Array state, i.e., the information of the availability of the task.

#### 4.2. Dynamics of the task allocation as a CPN

Typically, in these tree structures, every reference to a task is allocated on the leaves. Thus, every internal vertex allows to a working process taking a decision based on a function over the number of tasks in every leaf of the subtrees from each child of the current vertex. Hence, the decision function, which could be a simple deterministic algorithm, or a sophisticated statistical or stochastic procedure, is a key element for which a process traverses through the tree structure, in this kind of allocation techniques. The following definitions allow capturing these dynamics aspects of the task allocation based on a BDT.


Notice that, in Figure 2, the identification of a vertex v is given by vi,j where i ¼ depthðvÞ and j ¼ posðvÞ.

Definition 7. The w is a weight function defined from the set of edges E of a tree T to the set containing 0 and 1 as w : E ! f0, 1g, such that wðeÞ ¼ 0 if e connects a parent vertex to its left child, and 1 if it connects to its right child. Let consider a vertex v of a tree T denoted as vi,j such that i ¼ depthðvÞ and j ¼ posðvÞ. Then, for vi,x parent of left child vj, <sup>y</sup> and right child vj, <sup>z</sup>, with their respective edges ei,j,y and ei,j, <sup>z</sup>, with y < z, by definition of depth, then wðei,j,yÞ ¼ 0 and wðei,j, <sup>z</sup>Þ ¼ 1.

Notice that, the weight function w provides a zero for an edge connecting a father with its left child and a one for the edge connecting it with the right child. In the modeling methodology here introduced, it is assumed that there exists a data structure called Task Array, which may have a reference to a task, if it is initially inserted, and the task has not been taken yet. This Task Array is mapped to the location of the leaves of a tree of height H from left to right, i.e., the array has a capacity for referencing a maximum of The quaternion semi � widel tasks. Figure 1 shows a tree of height 3, with 8 leaves that are associated with the location in the Task Array, where references to task are T0, T1…, T7. Accordingly, the following functions are defined.


• The initial marking, assuming that the Task Array is initially full, is M<sup>0</sup> defined as follows:

Figure 4 shows the relation of the tree structure with the proposed CPN as stated by the previous definition. Observe that places of type O represent every level from the root to the leaves, forward, and backwards, i.e., a token in one of those place contains the information of pos, depth, and val. Places of type R represent the state of every subtree from every node at that level, i.e., the value of the data function for every node of the same level given by pos and depth. Place of

Typically, in these tree structures, every reference to a task is allocated on the leaves. Thus, every internal vertex allows to a working process taking a decision based on a function over the number of tasks in every leaf of the subtrees from each child of the current vertex. Hence, the decision function, which could be a simple deterministic algorithm, or a sophisticated statistical or stochastic procedure, is a key element for which a process traverses through the tree structure, in this kind of allocation techniques. The following definitions allow capturing

Definition 5. The pos is a position function defined from a set of vertexes V of a tree T to the set of natural numbers N as pos : V ! N, such that posðvÞ ¼ j, where j is the numeric position of vertex v with respect to all the vertexes in the same level and counting from left to right. Definition 6. The depth is a function defined from a set of vertexes V of a tree T to the set of natural numbers N as depth : V ! N, such that depthðvÞ ¼ #edges from the root node to v: Notice that, in Figure 2, the identification of a vertex v is given by vi,j where i ¼ depthðvÞ and

Definition 7. The w is a weight function defined from the set of edges E of a tree T to the set containing 0 and 1 as w : E ! f0, 1g, such that wðeÞ ¼ 0 if e connects a parent vertex to its left child, and 1 if it connects to its right child. Let consider a vertex v of a tree T denoted as vi,j such that i ¼ depthðvÞ and j ¼ posðvÞ. Then, for vi,x parent of left child vj, <sup>y</sup> and right child vj, <sup>z</sup>, with their respective edges ei,j,y and ei,j, <sup>z</sup>, with y < z, by definition of depth, then

Notice that, the weight function w provides a zero for an edge connecting a father with its left child and a one for the edge connecting it with the right child. In the modeling methodology

type A represents the Task Array state, i.e., the information of the availability of the task.

�2þ<sup>j</sup> <sup>¼</sup> 〈i, j, <sup>2</sup><sup>H</sup>�<sup>i</sup>

�<sup>1</sup>r2<sup>i</sup>

• M0½2ðH þ 1Þ þ H þ 1� ¼ <sup>∪</sup> <sup>i</sup>¼<sup>1</sup>nTai where ai ¼ 〈posi ¼ i, existi ¼ 1〉.

i∈ ½1, d�.

222 Computer Simulation

j ¼ posðvÞ.

wðei,j,yÞ ¼ 0 and wðei,j, <sup>z</sup>Þ ¼ 1.

• <sup>M</sup>0½jPOj þ <sup>i</sup>� ¼ <sup>∪</sup> <sup>i</sup>¼1<sup>H</sup> <sup>∪</sup> <sup>j</sup>¼02<sup>i</sup>

• M0½i� ¼ ∅ for i ∈½2, jPOj�.

4.2. Dynamics of the task allocation as a CPN

these dynamics aspects of the task allocation based on a BDT.

full, as mentioned.

• M0½1�¼fo1, …, odg, where oi ¼ depthi ¼ 0, posi ¼ 0, vali ¼ 0, successi ¼ 0 ¼ 0, 0, 0, 0 for

〉 assuming that the Task Array is initially

• (Induction): valðvi,jÞ ¼ val parentðvi,jÞ þ w eposðparentðvi,jÞ, depthðvi,jÞ,j �2<sup>H</sup>�<sup>i</sup> , for every pair ði, jÞ, where 1 ≤ i ≤ H, 0 ≤ j < 2<sup>i</sup> .

In Figure 2, for example, valðvÞ is shown in the left of every vertex v as valðvÞ=, e.g., valðv2,3Þ ¼ val parentðv2,3Þ þ w eposðparentðv2,3Þ, depthðv2,3Þ, <sup>3</sup> �23�<sup>2</sup> <sup>¼</sup> valðv1, <sup>1</sup>Þ þ <sup>w</sup> eposðv1, <sup>1</sup>Þ,2,<sup>3</sup> �21 <sup>¼</sup> val parentðv1,1Þ þ w eposðparentðv1,1Þ, depthðv1,1Þ, <sup>1</sup> �23�<sup>1</sup> <sup>þ</sup> <sup>w</sup> eposðv1,1Þ, <sup>2</sup>,<sup>3</sup> �2<sup>1</sup> <sup>¼</sup> valðv0,0Þþ¼ <sup>0</sup> þwðe0,1,1Þ�4 þ wðe1,2, <sup>3</sup>Þ�2 ¼ 0 þ 1�4 þ 1�2 ¼ 6.

Observe that function val generates the value of a node based on the value of his parent. Notice that, there are 8 leaves in level 3, and then, every node value is exactly one entry in the Task Array. For example, the value 4= at v3, <sup>5</sup> points to the entry T4. Now, consider the nodes in level 2, which are 4, exactly the half of those at level 3. Every vertex at this level, points to the first entry in a range of 2, i.e., v2, <sup>0</sup> has a value of 0, v2, <sup>1</sup> and has a value of 2. Thus, every vertex at level 1, points to the first entry of the two partitions of the Task Array, and therefore, v0, <sup>0</sup> points to the first entry of the whole of the Task Array.

For illustration purposes, suppose that there is one node with two threads and eight tasks to be executed. Accordingly, d ¼ 2 and nT ¼ 8. The CPN System ðN , M0Þ for modeling this problem is shown in Figure 4.

Notice that accordingly to the proposed method, a full binary tree is obtained when nT <sup>¼</sup> <sup>2</sup><sup>n</sup> for some n, since, in this case, H ¼ log2nT is exactly equal to log2nT. In the particular example, <sup>H</sup> <sup>¼</sup> 3, since nT <sup>¼</sup> 23 <sup>¼</sup> 8. Then, <sup>j</sup>POj ¼ <sup>2</sup>ð<sup>H</sup> <sup>þ</sup> <sup>1</sup>Þ ¼ 8, with PO ¼ fp1, p2, …, p8g; <sup>j</sup>PRj ¼ <sup>H</sup> <sup>¼</sup> 3, with PR ¼ fp9, p10, p11g; jPAj ¼ 1 with PA ¼ {p12}; jTj¼jPOj ¼ 8, with T ¼ {t1, t2, …, t8}; o ¼ {o1, <sup>o</sup>2} since <sup>d</sup> <sup>¼</sup> 2; <sup>r</sup> ¼ fr1, r2, …, r14g, since 2<sup>H</sup>þ<sup>1</sup> � <sup>2</sup> <sup>¼</sup> 23 � <sup>2</sup> <sup>¼</sup> 14, and <sup>a</sup> <sup>¼</sup> {a1, <sup>a</sup>2, …, <sup>a</sup>8}.

The following functions complement the CPN model:

• oi ¼ 〈depthi, posi, vali, successi〉 ¼ 〈0, 0, 0, 0〉, for i ∈{1, 2}. This marking represents the current state of every process, i.e., at the beginning, processes are at root node v0, <sup>0</sup> with a value 0.

• rj ¼ 〈depthj, posj, dataj〉, such that j ∈½1, jrj�, where: r<sup>1</sup> ¼ 〈1, 0, 4〉, r<sup>2</sup> ¼ 〈1, 1, 4〉, r<sup>3</sup> ¼ 〈2, 0, 2〉, r<sup>4</sup> ¼ 〈2, 1, 2〉, r<sup>5</sup> ¼ 〈2, 2, 2〉, r<sup>6</sup> ¼ 〈2, 3, 2〉, r<sup>7</sup> ¼ 〈3, 0, 1〉, r<sup>8</sup> ¼ 〈3, 1, 1〉, r<sup>9</sup> ¼ 〈3, 2, 1〉, r<sup>10</sup> ¼ 〈3, 3, 1〉, r<sup>11</sup> ¼ 〈3, 4, 1〉, r<sup>12</sup> ¼ 〈3, 5, 1〉, r<sup>13</sup> ¼ 〈3, 6, 1〉, r<sup>14</sup> ¼ 〈3, 7, 1〉. These markings provide the information about the data available for every subtree constructed from the node identified

a<sup>5</sup> ¼ 〈5, 1〉, a<sup>6</sup> ¼ 〈6, 1〉, a<sup>7</sup> ¼ 〈7, 1〉, a<sup>8</sup> ¼ 〈8, 1〉. These markings provide the information about the availability of the task located at pos, i.e., ty ¼ 1 if the task is available, 0 otherwise. The color mapping is CpðpkÞ ¼ o for k∈ ½1, 8�, CpðpkÞ ¼ r for k ∈½9, 11� and Cpðp12Þ ¼ a. That is, the places p1, …, p<sup>8</sup> accept tokens of type O, the places p9, p10, p<sup>11</sup> accept tokens of type R, and p<sup>12</sup> accept tokens of type A. The conditional color mapping is CtðtkÞ ¼ f <sup>1</sup> for k ∈½1, 3�, CtðtkÞ ¼ f <sup>2</sup> for k ∈½5, 7�, and CtðtkÞ ¼ f <sup>3</sup> for k ¼ 4 and k ¼ 9. The Pre- and Post matrices for the

Thus, Figure 4 represents a CPN that captures the structure of a BDT and the behavior of the working threads over the tree of height H ¼ log2nT ¼ log28 ¼ 3. They represents the places and transitions for the traveling from the root down to a leaf, and the reverse way from the leaf up to the root, on the left and right side of the CPN, respectively. Additionally, the central places, marked as "R," represent the "resources" in the system, i.e., the shared memory

Consider an initial token in place p<sup>1</sup> ¼ 〈0, 0, 0, 0〉 as described by the initial marking M0. Now, the binding for t<sup>1</sup> requires two tokens from p9, besides one initial token in p1, subject to

〉, such that l∈½1, nT�, where: a<sup>1</sup> ¼ 〈1, 1〉, a<sup>2</sup> ¼ 〈2, 1〉, a<sup>3</sup> ¼ 〈3, 1〉, a<sup>4</sup> ¼ 〈4, 1〉,

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

225

Success being 0 means that no task has been assigned to the current process.

as vdepth, pos.

• al ¼ 〈posl, tyl

registers and the tasks.

net of this example are shown in Figure 5.

Figure 5. Pre and postmatrices for CPN DS nT ¼ 8.

Figure 4. DS CPN with nT ¼ 8.


The color mapping is CpðpkÞ ¼ o for k∈ ½1, 8�, CpðpkÞ ¼ r for k ∈½9, 11� and Cpðp12Þ ¼ a. That is, the places p1, …, p<sup>8</sup> accept tokens of type O, the places p9, p10, p<sup>11</sup> accept tokens of type R, and p<sup>12</sup> accept tokens of type A. The conditional color mapping is CtðtkÞ ¼ f <sup>1</sup> for k ∈½1, 3�, CtðtkÞ ¼ f <sup>2</sup> for k ∈½5, 7�, and CtðtkÞ ¼ f <sup>3</sup> for k ¼ 4 and k ¼ 9. The Pre- and Post matrices for the net of this example are shown in Figure 5.

Thus, Figure 4 represents a CPN that captures the structure of a BDT and the behavior of the working threads over the tree of height H ¼ log2nT ¼ log28 ¼ 3. They represents the places and transitions for the traveling from the root down to a leaf, and the reverse way from the leaf up to the root, on the left and right side of the CPN, respectively. Additionally, the central places, marked as "R," represent the "resources" in the system, i.e., the shared memory registers and the tasks.

Consider an initial token in place p<sup>1</sup> ¼ 〈0, 0, 0, 0〉 as described by the initial marking M0. Now, the binding for t<sup>1</sup> requires two tokens from p9, besides one initial token in p1, subject to

Figure 5. Pre and postmatrices for CPN DS nT ¼ 8.

Figure 4. DS CPN with nT ¼ 8.

224 Computer Simulation

conditions given by f <sup>1</sup>: ðoi, rj, rkÞjdepthi þ 1 ¼ depthj ¼ depthk ∧ posi � 2 ¼ posj ∧ posi � 2 þ 1 ¼ posk. Since in this case oi ¼ o<sup>1</sup> ¼ 〈pos<sup>1</sup> ¼ 0, depht<sup>1</sup> ¼ 0, val<sup>1</sup> ¼ 0, success<sup>1</sup> ¼ 0〉 due to si1, then the required tokens from <sup>p</sup><sup>9</sup> are rj <sup>¼</sup> <sup>r</sup><sup>1</sup> <sup>¼</sup> 〈depth<sup>1</sup> <sup>¼</sup> <sup>0</sup>, pos<sup>1</sup> <sup>¼</sup> <sup>1</sup>, data<sup>1</sup> <sup>¼</sup> <sup>2</sup><sup>H</sup>�<sup>1</sup> <sup>¼</sup> <sup>4</sup>〉 and rk <sup>¼</sup> <sup>r</sup><sup>2</sup> <sup>¼</sup> 〈depth<sup>2</sup> <sup>¼</sup> <sup>1</sup>, pos<sup>2</sup> <sup>¼</sup> <sup>1</sup>, data<sup>2</sup> <sup>¼</sup> 23�<sup>1</sup> <sup>¼</sup> <sup>4</sup>〉, where data<sup>1</sup> <sup>¼</sup> data<sup>2</sup> <sup>¼</sup> 4, assuming that the Task Array is full of tasks to be attended. Thus, the output of t<sup>1</sup> due to so<sup>1</sup> is o<sup>1</sup> ¼ 〈depth<sup>1</sup> ¼ depth<sup>1</sup> þ 1, pos<sup>1</sup> <sup>¼</sup> pos<sup>1</sup> � <sup>2</sup> <sup>þ</sup> <sup>h</sup> <sup>¼</sup> <sup>1</sup>, val<sup>1</sup> <sup>þ</sup> <sup>h</sup> � 22 <sup>¼</sup> <sup>4</sup>, success<sup>1</sup> <sup>¼</sup> <sup>0</sup>〉 <sup>¼</sup> 〈1, <sup>1</sup>, <sup>4</sup>, <sup>0</sup>〉, assuming <sup>h</sup> <sup>¼</sup> 1.

as shown in (f), since the task reference represented by the marking < 6, 1 > is available as shown in (e). Notice that, this marking is updated to < 6, 0 > as shown in (f), as expected. In (g), the input state of the transitions subject to f <sup>2</sup> is represented with the marking < 3, 6, 6, 1 >. Then, it is binding by f <sup>2</sup> with the token < 3, 6, 1 >, as illustrated in (h). Thus, the input token < 3, 6, 6, 1 > is updated to < 2, 3, 6, 1 > by so2, as well as, < 3, 6, 1 > is updated to < 3, 6, 0 > by

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

227

One of the main advantages of modeling the task allocation problem by using CPN's is the possibility of varying the parameters during the simulation process, as well as the flexibility of the net structure in order to cope with a greater number of task to be attended. One of the key parameters for the simulation of the problem is the number of tokens of type O, which represents the number of processes or threads in a specific problem. The other important parameter is the decision function α, which is related to the spreading of processes or threads through the tree, by updating of depth and pos functions. It is clear, for example, that the distribution of the processes or threads at every level in the tree structure directly influences

The next section explains the simulation process that is possible to execute with the proposed modeling methodology and how these simulations can help explore different parameters in

This section shows how to simulate a CPN model as that obtained by the methodology detailed in the previous subsection. The simulation process allows to investigate the performance of the model under parametric variations. First, the general characteristics of the simulation framework are introduced. Then, the simulation and results of a CPN model, representing a problem with

In order to simplify the construction of the CPN model by using the herein proposed methodology, and thus, simulate its behavior, a CPN tool is used [43]. This environment allows editing, simulating and analyzing CPN models [36–42] and was successfully used for a variety

The model parameters that have to be considered are a number of initially available tasks ðnTasksÞ, the height of the tree with respect to the number of tasks ðHÞ, and the decision function ðαÞ. For this example, those values are nTasks ¼ 1024, H ¼ log21024 ¼ 10, and α ¼ ðr < dataj=ðdataj þ datakÞÞ for a given oi, rj, rk such that depthi þ 1 ¼ depthj ¼ depthk ∧ posi�2 ¼ posj ∧ posi�2 þ 1 ¼ posk and where r is a randomly generated number. These parameters are fixed during every simulation process. However, the number of threads is varying between

In the simulation framework, there are different measurements of parameters that can be obtained. For illustrative purposes, in this example, the attention is focused on study the impact on the performance of the task allocation procedure when the number of processes, or working

simulations, with <sup>d</sup> <sup>¼</sup> <sup>2</sup>n, <sup>1</sup> <sup>≤</sup> <sup>n</sup> <sup>≤</sup> nTasks, for a total of 10 simulations runs.

ro2, as shown in the section (i) of the figure.

the contention at the acquisition of the tasks.

order to optimize the task allocation problem.

5. Simulation of the CPN model

1024 tasks, are given.

threads, is increasing.

of applications areas [44–48].

The sketches of CPN in Figure 6 show the three main marking evolutions in the CPN subject to the binding functions f <sup>1</sup>, f <sup>2</sup>, f <sup>3</sup>. The marking evolution subject to f <sup>1</sup> from the current state is shown in (a). In (b), the tokens oi, rj, rk are taken by the depicted transition, since f <sup>1</sup>: ðoi, rj, rkÞjdepthi þ 1 ¼ depthj ¼ depthk ∧ posi�2 ¼ posj ∧posi�2 þ 1 ¼ posk, then oi ¼< 1, 1, 0, 0 > , rj ¼< 2, 2, 2 > , rk ¼< 2, 3, 2 > : In (c), the transition has fired, then o<sup>1</sup> ¼< 2, 3, 6, 0 >, since so1ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi <sup>þ</sup> <sup>1</sup>, posi � <sup>2</sup> <sup>þ</sup> hi, vali þ ðhiÞ � <sup>2</sup><sup>H</sup>�depthiþ<sup>1</sup> , successi〉 assuming that hi ¼ 1. The tokens rj, rk remain the same, since success ¼ 0. In (d), the marking evolution subject to f <sup>3</sup> from the current marking < 3, 6, 6, 0 > is updated to the marking < 3, 6, 6, 1 >

Figure 6. Binding and firing of a CPN transition.

as shown in (f), since the task reference represented by the marking < 6, 1 > is available as shown in (e). Notice that, this marking is updated to < 6, 0 > as shown in (f), as expected. In (g), the input state of the transitions subject to f <sup>2</sup> is represented with the marking < 3, 6, 6, 1 >. Then, it is binding by f <sup>2</sup> with the token < 3, 6, 1 >, as illustrated in (h). Thus, the input token < 3, 6, 6, 1 > is updated to < 2, 3, 6, 1 > by so2, as well as, < 3, 6, 1 > is updated to < 3, 6, 0 > by ro2, as shown in the section (i) of the figure.

One of the main advantages of modeling the task allocation problem by using CPN's is the possibility of varying the parameters during the simulation process, as well as the flexibility of the net structure in order to cope with a greater number of task to be attended. One of the key parameters for the simulation of the problem is the number of tokens of type O, which represents the number of processes or threads in a specific problem. The other important parameter is the decision function α, which is related to the spreading of processes or threads through the tree, by updating of depth and pos functions. It is clear, for example, that the distribution of the processes or threads at every level in the tree structure directly influences the contention at the acquisition of the tasks.

The next section explains the simulation process that is possible to execute with the proposed modeling methodology and how these simulations can help explore different parameters in order to optimize the task allocation problem.

## 5. Simulation of the CPN model

conditions given by f <sup>1</sup>: ðoi, rj, rkÞjdepthi þ 1 ¼ depthj ¼ depthk ∧ posi � 2 ¼ posj ∧ posi � 2 þ 1 ¼ posk. Since in this case oi ¼ o<sup>1</sup> ¼ 〈pos<sup>1</sup> ¼ 0, depht<sup>1</sup> ¼ 0, val<sup>1</sup> ¼ 0, success<sup>1</sup> ¼ 0〉 due to si1, then the required tokens from <sup>p</sup><sup>9</sup> are rj <sup>¼</sup> <sup>r</sup><sup>1</sup> <sup>¼</sup> 〈depth<sup>1</sup> <sup>¼</sup> <sup>0</sup>, pos<sup>1</sup> <sup>¼</sup> <sup>1</sup>, data<sup>1</sup> <sup>¼</sup> <sup>2</sup><sup>H</sup>�<sup>1</sup> <sup>¼</sup> <sup>4</sup>〉 and rk <sup>¼</sup> <sup>r</sup><sup>2</sup> <sup>¼</sup> 〈depth<sup>2</sup> <sup>¼</sup> <sup>1</sup>, pos<sup>2</sup> <sup>¼</sup> <sup>1</sup>, data<sup>2</sup> <sup>¼</sup> 23�<sup>1</sup> <sup>¼</sup> <sup>4</sup>〉, where data<sup>1</sup> <sup>¼</sup> data<sup>2</sup> <sup>¼</sup> 4, assuming that the Task Array is full of tasks to be attended. Thus, the output of t<sup>1</sup> due to so<sup>1</sup> is o<sup>1</sup> ¼ 〈depth<sup>1</sup> ¼ depth<sup>1</sup> þ 1,

The sketches of CPN in Figure 6 show the three main marking evolutions in the CPN subject to the binding functions f <sup>1</sup>, f <sup>2</sup>, f <sup>3</sup>. The marking evolution subject to f <sup>1</sup> from the current state is shown in (a). In (b), the tokens oi, rj, rk are taken by the depicted transition, since f <sup>1</sup>: ðoi, rj, rkÞjdepthi þ 1 ¼ depthj ¼ depthk ∧ posi�2 ¼ posj ∧posi�2 þ 1 ¼ posk, then oi ¼< 1, 1, 0, 0 > , rj ¼< 2, 2, 2 > , rk ¼< 2, 3, 2 > : In (c), the transition has fired, then o<sup>1</sup> ¼< 2, 3, 6, 0 >, since

ing that hi ¼ 1. The tokens rj, rk remain the same, since success ¼ 0. In (d), the marking evolution subject to f <sup>3</sup> from the current marking < 3, 6, 6, 0 > is updated to the marking < 3, 6, 6, 1 >

, successi〉 assum-

pos<sup>1</sup> <sup>¼</sup> pos<sup>1</sup> � <sup>2</sup> <sup>þ</sup> <sup>h</sup> <sup>¼</sup> <sup>1</sup>, val<sup>1</sup> <sup>þ</sup> <sup>h</sup> � 22 <sup>¼</sup> <sup>4</sup>, success<sup>1</sup> <sup>¼</sup> <sup>0</sup>〉 <sup>¼</sup> 〈1, <sup>1</sup>, <sup>4</sup>, <sup>0</sup>〉, assuming <sup>h</sup> <sup>¼</sup> 1.

226 Computer Simulation

so1ð〈depthi, posi, vali, successi〉Þ ¼ 〈depthi <sup>þ</sup> <sup>1</sup>, posi � <sup>2</sup> <sup>þ</sup> hi, vali þ ðhiÞ � <sup>2</sup><sup>H</sup>�depthiþ<sup>1</sup>

Figure 6. Binding and firing of a CPN transition.

This section shows how to simulate a CPN model as that obtained by the methodology detailed in the previous subsection. The simulation process allows to investigate the performance of the model under parametric variations. First, the general characteristics of the simulation framework are introduced. Then, the simulation and results of a CPN model, representing a problem with 1024 tasks, are given.

In order to simplify the construction of the CPN model by using the herein proposed methodology, and thus, simulate its behavior, a CPN tool is used [43]. This environment allows editing, simulating and analyzing CPN models [36–42] and was successfully used for a variety of applications areas [44–48].

The model parameters that have to be considered are a number of initially available tasks ðnTasksÞ, the height of the tree with respect to the number of tasks ðHÞ, and the decision function ðαÞ. For this example, those values are nTasks ¼ 1024, H ¼ log21024 ¼ 10, and α ¼ ðr < dataj=ðdataj þ datakÞÞ for a given oi, rj, rk such that depthi þ 1 ¼ depthj ¼ depthk ∧ posi�2 ¼ posj ∧ posi�2 þ 1 ¼ posk and where r is a randomly generated number. These parameters are fixed during every simulation process. However, the number of threads is varying between simulations, with <sup>d</sup> <sup>¼</sup> <sup>2</sup>n, <sup>1</sup> <sup>≤</sup> <sup>n</sup> <sup>≤</sup> nTasks, for a total of 10 simulations runs.

In the simulation framework, there are different measurements of parameters that can be obtained. For illustrative purposes, in this example, the attention is focused on study the impact on the performance of the task allocation procedure when the number of processes, or working threads, is increasing.

Table 1 summarizes some results obtained from the simulations performed on the CPN tools framework for a model with model described previously. The number of tasks was constant, while the number of processes was increased by a power of two, represented in the first column. The second column represents the average of events required per process to complete all the tasks. It provides a measure of the speed by increasing the number of processes. The third column represents the total number of readings to the task array. The forth column represents the total number of read collisions, i.e., when a process reads a task that was previously attended by other process. The fifth column represents the average of the reads to the task array executed by each process. The sixth column shows the average of the failed reads or collisions executed by each process. The seventh column shows the proficiency of the execution by the total amount of processes. Finally, the eighth column represents the average of total reads, including the failed ones, to the task array required per task.

The plot in Figure 7 shows the average of event, in the simulation framework, that each process required for completing all the tasks in the array. Notice that, as the number of processes increases exponentially, the number of events per processes decreases almost logarithmically.

The plot in Figure 8 shows the average of reads that each process has performed to the task array for the completion of all the tasks. As in the previous figure, notice that the number of reads decreases almost logarithmically as the number of processes increases exponentially.

The plot in Figure 9 shows the average of collisions, i.e., a process's read of a task that was previously attended by other process, as the number of processes increases. The figure shows that the maximum number of average collisions occurs when the number of processes is 32 in these experiments. The collisions, or failed reads, are highly influenced by the decision rule α.

The plot in Figure 10 shows the proficiency of attending all the tasks in the array as the number of processes increases. The proficiency increases almost exponentially as the number of


processes increases to 256. The ideal behavior is to double the proficiency as the number of processes also double. However, the collisions at reading the tasks undermine this proficiency,

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

229

as expected. Notice that, the plot is logarithmic.

Figure 8. Average of reads per process.

Figure 7. Average of events per process.

Table 1. Simulation results for nTasks ¼ 1024.

Figure 7. Average of events per process.

Table 1 summarizes some results obtained from the simulations performed on the CPN tools framework for a model with model described previously. The number of tasks was constant, while the number of processes was increased by a power of two, represented in the first column. The second column represents the average of events required per process to complete all the tasks. It provides a measure of the speed by increasing the number of processes. The third column represents the total number of readings to the task array. The forth column represents the total number of read collisions, i.e., when a process reads a task that was previously attended by other process. The fifth column represents the average of the reads to the task array executed by each process. The sixth column shows the average of the failed reads or collisions executed by each process. The seventh column shows the proficiency of the execution by the total amount of processes. Finally, the eighth column represents the average

The plot in Figure 7 shows the average of event, in the simulation framework, that each process required for completing all the tasks in the array. Notice that, as the number of processes increases exponentially, the number of events per processes decreases almost logarithmically.

The plot in Figure 8 shows the average of reads that each process has performed to the task array for the completion of all the tasks. As in the previous figure, notice that the number of reads decreases almost logarithmically as the number of processes increases exponentially.

The plot in Figure 9 shows the average of collisions, i.e., a process's read of a task that was previously attended by other process, as the number of processes increases. The figure shows that the maximum number of average collisions occurs when the number of processes is 32 in these experiments. The collisions, or failed reads, are highly influenced by the decision rule α. The plot in Figure 10 shows the proficiency of attending all the tasks in the array as the number of processes increases. The proficiency increases almost exponentially as the number of

Reads Collisions Reads per

 22528.00 1024.00 0.00 1024.00 0.00 1.00 1.00 11293.25 1027.00 3.00 513.50 1.50 1.99 1.00 5667.00 1030.67 6.67 257.67 1.67 3.97 1.01 2853.02 1038.00 14.00 129.75 1.75 7.89 1.01 1443.81 1051.33 27.33 65.71 1.71 15.58 1.03 743.34 1081.67 57.67 33.80 1.80 30.29 1.06 384.67 1117.33 93.33 17.46 1.46 58.65 1.09 202.77 1175.00 151.00 9.18 1.18 111.55 1.15 108.85 1265.67 241.67 4.94 0.94 207.12 1.24 55.39 1271.00 247.00 2.48 0.48 412.50 1.24 27.70 1285.50 261.50 1.26 0.26 815.70 1.26

process

Fails per process

Proficiency Read

per task

of total reads, including the failed ones, to the task array required per task.

Processes (d) Average

228 Computer Simulation

events

Table 1. Simulation results for nTasks ¼ 1024.

Figure 8. Average of reads per process.

processes increases to 256. The ideal behavior is to double the proficiency as the number of processes also double. However, the collisions at reading the tasks undermine this proficiency, as expected. Notice that, the plot is logarithmic.

Figure 9. Average of collisions per process.

6. Conclusions

Figure 11. Average of reads per task.

Author details

Mexico

construction and its respective simulation process.

\*Address all correspondence to: mildreth@iteso.mx

This work presented a framework, based on Colored Petri nets, for the modeling and simulation of task allocation problems, which arises in environments such as grid or cluster of computers. The framework allows the construction of complex problems in a compact way. Additionally, a simulation process of the constructed models permits the study of different key aspects of the allocation strategies. Some analysis that could be performed includes the impact of different decision rules of the processes for allocating the tasks, the effect of a greater number of processes in the contention of the task's acquisition or the ration of the increased speed to the attention of tasks by a greater number of processes, among others. Additionally, the methodology allows with ease the construction of structures for an incremental model

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

231

The proposed methodology allows with ease the extension to n ary tree structures as well as

tree structures with a greater number of tasks and their respective simulation process.

Mildreth Alcaraz-Mejia\*, Raul Campos-Rodriguez and Marco Caballero-Gutierrez

Electronics, Systems and Informatics Department, ITESO University, Tlaquepaque, Jalisco,

Figure 10. Average of proficiency.

The plot in Figure 11 shows the average of work that each process has to do for completing a task, measured as the number of reads. This number slightly increases as the number of processes increases due to the growth in the number of collisions.

Figure 11. Average of reads per task.

## 6. Conclusions

This work presented a framework, based on Colored Petri nets, for the modeling and simulation of task allocation problems, which arises in environments such as grid or cluster of computers. The framework allows the construction of complex problems in a compact way. Additionally, a simulation process of the constructed models permits the study of different key aspects of the allocation strategies. Some analysis that could be performed includes the impact of different decision rules of the processes for allocating the tasks, the effect of a greater number of processes in the contention of the task's acquisition or the ration of the increased speed to the attention of tasks by a greater number of processes, among others. Additionally, the methodology allows with ease the construction of structures for an incremental model construction and its respective simulation process.

The proposed methodology allows with ease the extension to n ary tree structures as well as tree structures with a greater number of tasks and their respective simulation process.

## Author details

The plot in Figure 11 shows the average of work that each process has to do for completing a task, measured as the number of reads. This number slightly increases as the number of

processes increases due to the growth in the number of collisions.

Figure 9. Average of collisions per process.

230 Computer Simulation

Figure 10. Average of proficiency.

Mildreth Alcaraz-Mejia\*, Raul Campos-Rodriguez and Marco Caballero-Gutierrez

\*Address all correspondence to: mildreth@iteso.mx

Electronics, Systems and Informatics Department, ITESO University, Tlaquepaque, Jalisco, Mexico

## References

[1] Lee, E. Y. S., and Tsuchiya, M., A task allocation model for distributed computing systems. IEEE Transactions on Computers, 1982, vol 100, no 1, pp. 41–47.

[14] Longo, F., Ghosh, R., Naik, V. K., and Trivedi, K. S., A scalable availability model for infrastructure-as-a-service cloud. In 2011 IEEE/IFIP 41st International Conference on

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

233

[15] Wei, B., Lin, C., and Kong, X., Dependability modeling and analysis for the virtual data center of cloud computing. In 2011 IEEE 13th International Conference on High Perfor-

[16] Li, X., Fan, Y., Sheng, Q. Z., Maamar, Z., and Zhu, H., A petri net approach to analyzing behavioral compatibility and similarity of web services. IEEE Transactions on Systems, Man and Cybernetics, Part A: Systems and Humans, 2011, vol. 41, no 3, pp. 510–521. [17] Muñoz, D. M., Correcher, A., García, E., and Morant, F., Stochastic DES fault diagnosis with coloured interpreted Petri nets. Mathematical Problems in Engineering, vol. 2015,

[18] Ghainani, A. T., Mohd Zin, A. A., and Ismail, N. A. M., Fuzzy timing Petri net for fault diagnosis in power system. Mathematical Problems in Engineering, vol. 2012, Article ID

[19] Alcaraz-Mejia, M., Lopez-Mellado, E., Ramirez-Treviño, A., Rivera-Rangel, I., Petri net based fault diagnosis of discrete event systems, In Proceedings of the IEEE International

[20] Latorre, J. I., Jiménez, E., and Pérez, M., The optimization problem based on alternatives aggregation Petri nets as models for industrial discrete event systems, Simulation, March

[21] Latorre, J. I., and Jiménez, E., Simulation-based optimization of discrete event systems with alternative structural configurations using distributed computation and the Petri net

[22] Alcaraz-Mejia, M., and Lopez-Mellado, E., Petri net model reconfiguration of discrete manufacturing systems. In Alexandre Dolgui Gerard Morel Carlos Pereira Editors: Information Control Problems in Manufacturing. 1st edition, Editorial Elsevier Science, 2006. pp. 517–552. eBook ISBN: 9780080478487, ISBN: 9780080446547, Page Count:

[23] Alcaraz-Mejia, M., and Lopez-Mellado, E., Fault recovery of manufacturing systems based on controller reconfiguration. In 2006 IEEE/SMC International Conference on Sys-

[24] Alcaraz-Mejia, M., Lopez-Mellado, E., and Ramirez-Treviño, A., A redundancy based method for Petri net model reconfiguration. In IEEE International Conference on Sys-

[25] Chen, S. J., Zhan, T. S., Huang, C. H., Chen, J. L., and Lin, C. H. (2015). Nontechnical loss and outage detection using fractional-order self-synchronization error-based fuzzy petri nets in micro-distribution systems. IEEE Transactions on smart grid, 6(1), 411–420.

Conference on Systems, Man and Cybernetics. pp. 4730–4735, October 2003.

paradigm, Simulation, November 2013; vol. 89, no 11: pp. 1310–1334.

mance Computing and Communications (HPCC), IEEE, 2011. pp. 784–789.

Dependable Systems and Networks (DSN), IEEE 2011. pp. 335–346.

Article ID 303107, 13 pages, 2015. doi:10.1155/2015/303107.

717195, 12 pages, 2012. doi:10.1155/2012/717195.

tem of Systems Engineering, IEEE, 2006. p. 6.

tems, Man and Cybernetics, 2007. IEEE, 2007. pp. 1382–1387.

2013; vol. 89, no 3: pp. 346–361.

2480.


[14] Longo, F., Ghosh, R., Naik, V. K., and Trivedi, K. S., A scalable availability model for infrastructure-as-a-service cloud. In 2011 IEEE/IFIP 41st International Conference on Dependable Systems and Networks (DSN), IEEE 2011. pp. 335–346.

References

232 Computer Simulation

[1] Lee, E. Y. S., and Tsuchiya, M., A task allocation model for distributed computing sys-

[2] Han, Y., Jiang, C., and Luo, X., Resource scheduling model for grid computing based on sharing synthesis of Petri net, in Proceedings of the 9th International Conference on Computer Supported Cooperative Work in Design, Coventry, UK, 2005, pp. 367–372. [3] Han, Y., Jiang, C., and Luo, X., Modelling and performance analysis of grid task scheduling based on composition and reduction of Petri nets, in Proceedings of the 5th International Conference on Grid and Cooperative Computing, Changsha, China, 2006, pp. 331–334. [4] Zhao, X., Wang, B., and Xu, L., Grid application scheduling model based on Petri net with changeable structure, in Proceeding of 6th International Conference on Grid and Cooper-

[5] Han, Y., Jiang, C., and Luo, X., Resource scheduling scheme for grid computing and its Petri net model and analysis, Parallel and Distributed Processing and Applications, Lecture Notes in Computer Science, vol. 3759, G. Chen et al., editors, Heidelberg:

[6] Hu, Z. G., Hu, R., Gui, W. H., Chen, J. E., and Chen, S. Q., General scheduling framework in computational grid based on Petri net, Journal of Central South University of Technol-

[7] Shojafar, M., Barzegar, S., and Meybodi, M. R., A new method on resource scheduling in grid systems based on hierarchical stochastic Petri net, in Proceedings of the 3rd International Conference on Computer and Electrical Engineering, Chengdu, China, 2010, pp. 175–180.

[8] Buyya, R., and Murshed, M., GridSim: a toolkit for the modelling and simulation of distributed resource management and scheduling for grid computing, Concurrency and

[9] Zhovtobryukh, D., A Petri net-based approach for automated goal-driven web service

[10] Haggarty, O. J., Knottenbelt, W. J., and Bradley, J. T., Distributed response time analysis of GSPN models with MapReduce, Simulation, August 2009; vol. 85, 8: pp. 497–509.

[11] Narciso, M., Piera, M. A., and Guasch, A., A methodology for solving logistic optimization problems through simulation, Simulation, May/June 2010; vol. 86, 5–6: pp. 369–389.

[12] Xiong, Z. G., Zhai, Z. L., Zhang, X. M., and Xia, X. W., Grid workflow service composition based on colored petri net. JDCTA: International Journal of Digital Content Technology

[13] Camilli, M., Petri nets state space analysis in the cloud. In 2012 34th International Confer-

Computation: Practice and Experience, vol. 14, no. 13–15, pp. 1175–1220, 2002.

composition, Simulation, January 2007; vol. 83, 1: pp. 33–63.

and its Applications, 2011, vol. 5, no 5, pp. 125–131.

ence on Software Engineering (ICSE). pp. 1638–1640.

tems. IEEE Transactions on Computers, 1982, vol 100, no 1, pp. 41–47.

ative Computing, Los Alamitos, CA, 2007, pp. 733–736.

Springer, pp. 530–539, 2005.

ogy, vol. 12, no. 1, pp. 232–237, 2005.


[26] Muñoz, D. M., Correcher, A., García, E., and Morant, F., Identification of stochastic timed discrete event systems with st-IPN. Mathematical Problems in Engineering, vol. 2014, Article ID 835312, 21 pages, 2014. doi:10.1155/2014/835312.

[39] Wells, L., Performance analysis using CPN tools. In Proceedings of the 1st international conference on Performance evaluation methodologies and tools. ACM, 2006. p. 59. [40] Westergaard, M., CPN Tools 4: Multi-formalism and extensibility. In Application and Theory of Petri Nets and Concurrency. Springer Berlin Heidelberg, 2013. pp. 400–409. [41] Westergaard, M. and Slaats, T., CPN Tools 4: a process modeling tool combining declarative and imperative paradigms. Automatic Control and Computer Sciences, 2013, vol.

Modeling and Simulation of Task Allocation with Colored Petri Nets

http://dx.doi.org/10.5772/67950

235

[42] Dworzański, L. W., and Lomazova, Irina A., CPN Tools-assisted simulation and verification of nested Petri nets. Automatic Control and Computer Sciences, 2013, vol. 47, no 7, pp. 393–402.

[43] Cpntools.org, (2015). CPN Tools Homepage. [online] Available at: http://cpntools.org/

[44] Machado, R. J., Lassen, K. B., Oliveira, S., Couto, M., and Pinto, P., Requirements validation: execution of UML models with CPN tools. International Journal on Software Tools

[45] Vanderfeesten, I., van der Aalst, W., and Reijers, H. A., Modelling a product based workflow system in CPN tools. In Proceedings of the Sixth Workshop on the Practical

[46] Störrle, H., Semantics and verification of data flow in UML 2.0 activities. Electronic Notes

[47] Yi, X., and Kochut, K. J., A cp-nets-based design and verification framework for web services composition. In Web Services, 2004. Proceedings. IEEE International Conference

[48] Rozinat, A., Wynn, M. T., van der Aalst, W. M., terHofstede, A. H., and Fidge, C. J., Workflow simulation for operational decision support. Data and Knowledge Engineer-

Use of Coloured Petri Nets and CPN Tools (CPN 2005). 2005. pp. 99–118.

in Theoretical Computer Science, 2005, vol. 127, no 4, pp. 35–52.

for Technology Transfer, 2007, vol. 9, no 3–4, pp. 353–369.

47, no 7, pp. 393–402.

[Accessed 24 Jul. 2015].

on. IEEE, 2004. pp. 756–760.

ing, 2009, vol. 68, no 9, pp. 834–850.


[39] Wells, L., Performance analysis using CPN tools. In Proceedings of the 1st international conference on Performance evaluation methodologies and tools. ACM, 2006. p. 59.

[26] Muñoz, D. M., Correcher, A., García, E., and Morant, F., Identification of stochastic timed discrete event systems with st-IPN. Mathematical Problems in Engineering, vol. 2014,

[27] Alistarh, D., Bender, M., Gilbert, S., and Guerraoui, R., How to allocate tasks asynchronously. 53rd Annual IEEE Symposium on Foundations of Computer Science, FOCS 2012,

[28] Jevtić, A., Gutiérrez, A., Andina, D., and Jamshidi, M., Distributed bees algorithm for task allocation in swarm of robots. IEEE Systems Journal, 2012, vol. 6, no 2, pp. 296–304.

[29] Delle Fave, F. M., Rogers, A., Xu, Z., Sukkarieh, S., and Jennings, N. R., Deploying the max-sum algorithm for decentralised coordination and task allocation of unmanned aerial vehicles for live aerial imagery collection. In 2012 IEEE International Conference

[30] Johnson, L., Choi, H. L., Ponda, S., and How, J. P., Allowing non-submodular score functions in distributed task allocation. In 2012 IEEE 51st Annual Conference on Decision

[31] Zhao-Pin, S. U., Jiang, J. G., Liang, C. Y., and Zhang, G. F., A distributed algorithm for parallel multi-task allocation based on profit sharing learning. Acta Automatica Sinica,

[32] Macarthur, K. S., Stranders, R., Ramchurn, S. D., and Jennings, N. R., A distributed anytime algorithm for dynamic task allocation in multi-agent systems. In Proceedings of the 25th Association of the Advancement on Artificial Intelligence Conference (AAAI 2011), pp.

[33] Alistarh, D., Aspnes, J., Bender, M. A., Gelashvili, R., and Gilbert, S. Dynamic task allocation in asynchronous shared memory. In Proceedings of the Twenty-Fifth Annual

[34] Diestel, R., Graph Theory. Springer-Verlag, Heidelberg. Berlin. Graduate Texts in Mathe-

[35] Girault, C. and Valk, R., Petri Nets for Systems Engineering: A Guide to Modelling, Verification, and Applications, July 30, 2001, Springer-Verlag Berlin Heidelberg, Newyork, 2003.

[36] Jensen, K., Coloured Petri nets: basic concepts, analysis methods and practical use.

[37] Jensen, K., Kristensen, L. M., and Wells, L., Coloured Petri Nets and CPN Tools for modelling and validation of concurrent systems. International Journal on Software Tools

[38] Ratzer A.V. et al. (2003) CPN Tools for Editing, Simulating, and Analysing Coloured Petri Nets. In: van der Aalst W.M.P., Best E. (eds) Applications and Theory of Petri Nets 2003. ICATPN 2003. Lecture Notes in Computer Science, vol. 2679. Springer, Berlin, Heidelberg.

ACM-SIAM Symposium on Discrete Algorithms. SIAM, 2014. pp. 416–435.

ISBN: 3642074472 9783642074479. Series ISSN: 1431–2654. vol. 2, no 1.

for Technology Transfer, 2007, vol. 9, no 3–4, pp. 213–254.

Article ID 835312, 21 pages, 2014. doi:10.1155/2014/835312.

New Brunswick, NJ, USA, October 20-23, 2012, pp. 331–340, 2012.

on Robotics and Automation (ICRA). IEEE, 2012, pp. 469–476.

and Control (CDC). IEEE, 2012, pp. 4702–4708.

2011, vol. 37, no 7, pp. 865–872.

234 Computer Simulation

356–362, 2011. San Francisco, USA.

matics, vol. 173, 3rd edition, 2005.

Springer Science and Business Media, 2013.


**Chapter 11**

**Multi-Criteria Decision-Making in the Implementation**

The consideration of renewable energy sources as sources for the production of electricity, demands an approach that would enable an analysis which comprehends various factors and stakeholders. The *P*reference *R*anking *O*rganization *METH*od for *E*nrichment *E*valuations (PROMETHEE), as a mathematical model for multi-criteria decision-making, is one of the ideal methods used when it is necessary to rank scenarios according to specific criteria, depending on whom the ranking is applied. This chapter presents various scenarios whose ranking is done according to defined criteria and weight coefficients for each of the stakeholders. This model recognized and accepted according to the theory of decision-making could be used as a tool for so-called stakeholder value approach. **Keywords:** renewable energy sources, PROMETHEE, the production of electricity, stakeholder value, multi-criteria decision-making proces, National Renewable Energy Action Plan (NAPOIE), mini hydros, biomass, wind, solar, geothermal energy

The basis for this chapter was document which established the goals in usage of renewable energy sources until 2020 (National Renewable Energy Action Plan of the Republic of Serbia further on NAPOIE) [2], as well as the manner in which they are to be achieved. In addition,

'According to article 20 of the Treaty Establishing Energy Community (further on: UOEnZ), the Republic of Serbia accepted the obligation to apply European Directives in the field of renewable energy sources (further on: OIE) —Directive 2001/77/EC on the promotion of the use of energy from renewable sources and Directive 2003/30/EC on the promotion of the use

> © 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

it has the goal to enhance investments in the field of renewable energy sources.

**of Renewable Energy Sources**

Additional information is available at the end of the chapter

Dejan Jovanovic and Ivan Pribicevic

http://dx.doi.org/10.5772/67734

**Abstract**

**1. Introduction**

## **Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources**

Dejan Jovanovic and Ivan Pribicevic

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/67734

#### **Abstract**

The consideration of renewable energy sources as sources for the production of electricity, demands an approach that would enable an analysis which comprehends various factors and stakeholders. The *P*reference *R*anking *O*rganization *METH*od for *E*nrichment *E*valuations (PROMETHEE), as a mathematical model for multi-criteria decision-making, is one of the ideal methods used when it is necessary to rank scenarios according to specific criteria, depending on whom the ranking is applied. This chapter presents various scenarios whose ranking is done according to defined criteria and weight coefficients for each of the stakeholders. This model recognized and accepted according to the theory of decision-making could be used as a tool for so-called stakeholder value approach.

**Keywords:** renewable energy sources, PROMETHEE, the production of electricity, stakeholder value, multi-criteria decision-making proces, National Renewable Energy Action Plan (NAPOIE), mini hydros, biomass, wind, solar, geothermal energy

## **1. Introduction**

The basis for this chapter was document which established the goals in usage of renewable energy sources until 2020 (National Renewable Energy Action Plan of the Republic of Serbia further on NAPOIE) [2], as well as the manner in which they are to be achieved. In addition, it has the goal to enhance investments in the field of renewable energy sources.

'According to article 20 of the Treaty Establishing Energy Community (further on: UOEnZ), the Republic of Serbia accepted the obligation to apply European Directives in the field of renewable energy sources (further on: OIE) —Directive 2001/77/EC on the promotion of the use of energy from renewable sources and Directive 2003/30/EC on the promotion of the use

of biofuels and other renewable fuels for transport. Those Directives were gradually replaced since 2009, and in January 2012 abolished by the new Directive 2009/28/EC of the European Parliament and the Council of the 23rd of April 2009 on the promotion of the use of energy from renewable sources, amending and subsequently repealing Directives 2001/77/EC and 2003/30/EC CELEX No. 32009L0028'.1 \*\*

With the adoption of the 'Law of Ratification…',<sup>2</sup> [3] the Republic of Serbia internationally committed to create NAPOIE [2].

Data given in **Tables 1** and **2** were used as input data for this chapter.

Types of renewable energy sources taken in consideration in this chapter are as follows:


The National Renewable Energy Action Plan of the Republic of Serbia (NAPOIE) defined target values, that is, the amount of GWh expected to be produced from every renewable energy source and to be delivered in the system. The defined goal is 2252 GWh obtained from following renewable energy sources: mini hydros, biomass, solar, wind and geothermal energy (**Table 3**).

**Type of renewable energy sources (MW) (GWh) Specific investment** 

**Type of renewable energy sources (MW) Estimated work** 

Biomass: power plants with combined

Biomass: power plants with combined

Biogas (manure): power plants with

production

1

combined production

Biogas (manure): power plants with

production

1

combined production

HE (over 10 MW) 250 1108 1819 454.8 MHE (up to 10 MW) 188 592 2795 525.5 Plants powered by wind energy 500 1000 1417 708.5 Plants powered by solar energy 10 13 2500 25.0

**Table 1.** The production of electricity from renewable energy sources from new plants in 20201

Geothermal energy 1 7 4115 4.1 Waste 3 18 4147 12.4 Landfill gas 10 50 2000 20.0 Total planned capacity 1092 3653 – 2322.6

**costs\* (€/kW)**

 Full name "Law on Ratification of the Treaty Establishing Energy Community between the European Community and the Republic of Albania, Republic of Bulgaria, Bosnia and Herzegovina, Republic of Croatia, Former Yugoslav Republic of Macedonia, Republic of Montenegro, Romania, Republic of Serbia and United Nations Interim Administration Mission on Kosovo in compliance with the Resolution 1244 of the UN Security Council" ("Službeni glasnik RS", no. 62/06 [2].

**hours (h)**

100 6400 640 55 17.5

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

30 7500 225 19 6.2

HE (over 10 MW) 250 4430 1108 95 30.3 MHE (up to 10 MW) 188 3150 592 51 16.2 Wind energy 500 2000 1000 86 27.4 Solar energy 10 1300 13 1 0.4

Geothermal energy 1 7000 7 1 0.2 Waste 3 6000 18 2 0.5 Landfill gas 10 5000 50 4 1.4 Total planned capacity 1092 – 3653 314 100.0

100 640 4522 452.2

30 225 4006 120.2

 Full name "Law on Ratification of the Treaty Establishing Energy Community between the European Community and the Republic of Albania, Republic of Bulgaria, Bosnia and Herzegovina, Republic of Croatia, Former Yugoslav Republic of Macedonia, Republic of Montenegro, Romania, Republic of Serbia and United Nations Interim Administration Mission on Kosovo in compliance with the Resolution 1244 of the UN Security Council" ("Službeni glasnik RS", no. 62/06 [2].

**Table 2.** Estimated finances for each of the technologies using renewable energy sources in the production of electricity needed to complete the planned share in energy production from new capacities until 2020 in electric energy sector1

**Price according to planned installed capacity until 2020 (millions €)**

.

**(GWh) (ktoe) Participation (%)**

http://dx.doi.org/10.5772/67734

239

.

The *goal* is to verify the ranking sequence of renewable energy sources if only one of the listed renewable energy sources would be delivering the total expected amount of GWh into the system and to rank scenarios according to stakeholders3 , on the basis of previously defined criteria and calculated weight coefficients, and also to establish whether the sequence of renewable energy sources is identical for all stakeholders.

On the basis of ranking achieved this way, we may determine which type of renewable energy source is the priority, depending on the stakeholder, and also whether the participation of all listed types is justified.

A multi-criteria analysis will provide a clearly established sequence of renewable energy sources for the stakeholders, and according to clearly established criteria. This sequence is important for the establishing of priorities.

<sup>1</sup> Taken from introduction of document NAPOIE, Ministarstvo energetike, razvoja i zaštite životne sredine, strana 18, Beograd 2013 [2].

<sup>2</sup> Full name "Law on Ratification of the Treaty Establishing Energy Community between the European Community and the Republic of Albania, Republic of Bulgaria, Bosnia and Herzegovina, Republic of Croatia, Former Yugoslav Republic of Macedonia, Republic of Montenegro, Romania, Republic of Serbia and United Nations Interim Administration Mission on Kosovo in compliance with the Resolution 1244 of the UN Security Council" ("Službeni glasnik RS", no. 62/06). 3 R. Edward Freeman. The stakeholder theory is a theory of organizational management and business ethics that addresses morals and values in managing an organization [1].


1 Full name "Law on Ratification of the Treaty Establishing Energy Community between the European Community and the Republic of Albania, Republic of Bulgaria, Bosnia and Herzegovina, Republic of Croatia, Former Yugoslav Republic of Macedonia, Republic of Montenegro, Romania, Republic of Serbia and United Nations Interim Administration Mission on Kosovo in compliance with the Resolution 1244 of the UN Security Council" ("Službeni glasnik RS", no. 62/06 [2].

**Table 1.** The production of electricity from renewable energy sources from new plants in 20201 .

of biofuels and other renewable fuels for transport. Those Directives were gradually replaced since 2009, and in January 2012 abolished by the new Directive 2009/28/EC of the European Parliament and the Council of the 23rd of April 2009 on the promotion of the use of energy from renewable sources, amending and subsequently repealing Directives 2001/77/EC and

Types of renewable energy sources taken in consideration in this chapter are as follows:

The National Renewable Energy Action Plan of the Republic of Serbia (NAPOIE) defined target values, that is, the amount of GWh expected to be produced from every renewable energy source and to be delivered in the system. The defined goal is 2252 GWh obtained from following renewable energy sources: mini hydros, biomass, solar, wind and geothermal energy

The *goal* is to verify the ranking sequence of renewable energy sources if only one of the listed renewable energy sources would be delivering the total expected amount of GWh into the

criteria and calculated weight coefficients, and also to establish whether the sequence of

On the basis of ranking achieved this way, we may determine which type of renewable energy source is the priority, depending on the stakeholder, and also whether the participation of all

A multi-criteria analysis will provide a clearly established sequence of renewable energy sources for the stakeholders, and according to clearly established criteria. This sequence is

Taken from introduction of document NAPOIE, Ministarstvo energetike, razvoja i zaštite životne sredine, strana 18,

Full name "Law on Ratification of the Treaty Establishing Energy Community between the European Community and the Republic of Albania, Republic of Bulgaria, Bosnia and Herzegovina, Republic of Croatia, Former Yugoslav Republic of Macedonia, Republic of Montenegro, Romania, Republic of Serbia and United Nations Interim Administration Mission on Kosovo in compliance with the Resolution 1244 of the UN Security Council" ("Službeni glasnik RS", no. 62/06).

R. Edward Freeman. The stakeholder theory is a theory of organizational management and business ethics that ad-

[3] the Republic of Serbia internationally

, on the basis of previously defined

\*\*

Data given in **Tables 1** and **2** were used as input data for this chapter.

2003/30/EC CELEX No. 32009L0028'.1

committed to create NAPOIE [2].

• Mini hydros (up to 10 MW).

• Wind energy. • Solar energy.

• Geothermal energy.

listed types is justified.

• Biomass.

238 Computer Simulation

(**Table 3**).

1

2

3

Beograd 2013 [2].

With the adoption of the 'Law of Ratification…',<sup>2</sup>

system and to rank scenarios according to stakeholders3

renewable energy sources is identical for all stakeholders.

important for the establishing of priorities.

dresses morals and values in managing an organization [1].


1 Full name "Law on Ratification of the Treaty Establishing Energy Community between the European Community and the Republic of Albania, Republic of Bulgaria, Bosnia and Herzegovina, Republic of Croatia, Former Yugoslav Republic of Macedonia, Republic of Montenegro, Romania, Republic of Serbia and United Nations Interim Administration Mission on Kosovo in compliance with the Resolution 1244 of the UN Security Council" ("Službeni glasnik RS", no. 62/06 [2].

**Table 2.** Estimated finances for each of the technologies using renewable energy sources in the production of electricity needed to complete the planned share in energy production from new capacities until 2020 in electric energy sector1 .


**Table 3.** Available potentials [4].

For solving this type of problems, one of the mathematical models that can be used is the one developed by Jean-Pierre Brans in 1982, for a multi-criteria decision-making in a group of alternatives described with several attributes.

## **2. Theoretical overview of the PROMETHEE**

The *P*reference *R*anking *O*rganization *METH*od for *E*nrichment *E*valuations (PROMETHEE)4 is part of a group of methods for multi-criteria decision-making within a group of alternatives described with several attributes, used as criteria. This method enables a comprehensive structuring of quality and quantity criteria of different importance into a relation of partial organization in a unique result (PROMETHEE II), on the basis of which alternatives can be ranked in an absolute manner.

We will consider a multi-criteria problem:

$$\text{Max}\{\left(k\_1(a), \dots, k\_l(a)\right) \mid a \in A\},\tag{1}$$

between two alternatives (action and activity), according to every criterion, on the basis of the difference (differentiation) of criteria values of alternatives which are being compared.

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

http://dx.doi.org/10.5772/67734

241

PROMETHEE II can be a tool for 'Management philosophy that regards maximization of the interests of its all stakeholders (customers, employees, shareholders and the community) as

The preference relation obtained this way is used so that input and output flows are calculated for each alternative, in graphs or tables. On the basis of these flows, the decision maker can apply partial ranking (PROMETHEE I) or absolute ranking (PROMETHEE II) in

Let *k* be a real function used to express one of the attributes used as a criterion for comparing

*k* : *A* → **R** (2)

Let us assume that this is a usefulness criterion, that is, that alternatives (scenarios/models)

For every alternative *a* d*A, k*(*a*) a criterion value is calculated according to criterion *k*. When two alternatives *a, b* d*A* are being compared, the result of that comparison is expressed as

 *P* : *A* × *A* → [0, 1] (3) the intensity of preference for alternative *a* in relation to alternative *b* is expressed, with the

Preference function that is added to a given criterion is the difference function of criteria value

 *P*(*a*, *b* ) = *P*(*k*(*a* ) −*k*(*b* ) )= *P*(*d* ) (4) *P*(*d*) is a non-decreasing function that assumes value zero for negative difference values *d* = *k* (*a*) − *k* (*b*) , *if the functions should be maximized, that is, P* (*a, b*) = *P* (*k* (*a*) − *k* (*b*)), that is, *d* = −(*k*

*P* (*a, b*) = 0 marks indifference between *a* and b, that is, there is no preference of *a* over *b*,

are compared according to this criterion on the basis of the principle 'bigger is better'.

In this chapter, the absolute ranking method PROMETHEE II was used.

its highest objective'.<sup>5</sup>

the group of alternatives.

alternatives:

a preference.

5

With preference function *P*

following interpretation:

*P* (*a, b*) *≈* 0 marks weak preference of *a* over b,

*P* (*a, b*) *≈* 1 marks strong preference of *a* over b,

*P* (*a, b*) = 1 marks strict preference of *a* over *b*.

(*a*) − *k* (*b*)) if the criterion is minimized (**Table 4**).

http://www.businessdictionary.com/definition/stakeholder-value-approach.html.

of alternatives, and it can be written as

**2.1. PROMETHEE preference relation**

where *A* is a finite group of activities and *ki* = 1,…, *k* are *usefulness criteria* which should be maximized or fulfilled according to the principle 'bigger is better' (this supposition enables a more simple presentation of the method—in cases when some of the criteria are *price criteria*, they can be transformed into usefulness criteria, or we can adjust the proceeding to those criteria as well).

The application of the PROMETHEE is characterized by two steps:


In the first step, a complex preference relation is formed (in order to stress the fact that this relation is based on the consideration of more criteria, this relation is called *outranking relation*), based on the generalization of the notion of the criteria. A *preference index* is then defined and a complex preference relation is obtained, which is shown in a graph representation. The essence of this step is that the decision maker (stakeholder) must express his preference

<sup>4</sup> Theoretical overview of the PROMEHTEE method is described in brief according to the "*Odlučivanje"*, Milutin Čupić, Milija Suknović, Fakultet organizacionih nauka, Beograd 2010. All general theoretical formulas, functions and graphs are taken from Ref. [5].

between two alternatives (action and activity), according to every criterion, on the basis of the difference (differentiation) of criteria values of alternatives which are being compared.

PROMETHEE II can be a tool for 'Management philosophy that regards maximization of the interests of its all stakeholders (customers, employees, shareholders and the community) as its highest objective'.<sup>5</sup>

The preference relation obtained this way is used so that input and output flows are calculated for each alternative, in graphs or tables. On the basis of these flows, the decision maker can apply partial ranking (PROMETHEE I) or absolute ranking (PROMETHEE II) in the group of alternatives.

In this chapter, the absolute ranking method PROMETHEE II was used.

### **2.1. PROMETHEE preference relation**

For solving this type of problems, one of the mathematical models that can be used is the one developed by Jean-Pierre Brans in 1982, for a multi-criteria decision-making in a group of

The *P*reference *R*anking *O*rganization *METH*od for *E*nrichment *E*valuations (PROMETHEE)4 is part of a group of methods for multi-criteria decision-making within a group of alternatives described with several attributes, used as criteria. This method enables a comprehensive structuring of quality and quantity criteria of different importance into a relation of partial organization in a unique result (PROMETHEE II), on the basis of which alternatives can be ranked in

(*a* ) ,…,*kk*

The application of the PROMETHEE is characterized by two steps:

(2) using this relation to find an answer to the problem (1.1).

(1) constructing a preference relation within a group of alternatives *A*,

mized or fulfilled according to the principle 'bigger is better' (this supposition enables a more simple presentation of the method—in cases when some of the criteria are *price criteria*, they can be transformed into usefulness criteria, or we can adjust the proceeding to those criteria as well).

In the first step, a complex preference relation is formed (in order to stress the fact that this relation is based on the consideration of more criteria, this relation is called *outranking relation*), based on the generalization of the notion of the criteria. A *preference index* is then defined and a complex preference relation is obtained, which is shown in a graph representation. The essence of this step is that the decision maker (stakeholder) must express his preference

Theoretical overview of the PROMEHTEE method is described in brief according to the "*Odlučivanje"*, Milutin Čupić, Milija Suknović, Fakultet organizacionih nauka, Beograd 2010. All general theoretical formulas, functions and graphs

(*a* ))|*a* ∈ *A*}, (1)

= 1,…, *k* are *usefulness criteria* which should be maxi-

alternatives described with several attributes.

**Renewable energy type Mtoe** Hydro 0.80 Solar 0.60 Biomass 2.25 Wind 0.20 Geothermal energy 0.20

We will consider a multi-criteria problem:

where *A* is a finite group of activities and *ki*

Max{(*k*<sup>1</sup>

an absolute manner.

**Table 3.** Available potentials [4].

240 Computer Simulation

4

are taken from Ref. [5].

**2. Theoretical overview of the PROMETHEE**

Let *k* be a real function used to express one of the attributes used as a criterion for comparing alternatives:

$$k: A \to \mathbf{R} \tag{2}$$

Let us assume that this is a usefulness criterion, that is, that alternatives (scenarios/models) are compared according to this criterion on the basis of the principle 'bigger is better'.

For every alternative *a* d*A, k*(*a*) a criterion value is calculated according to criterion *k*. When two alternatives *a, b* d*A* are being compared, the result of that comparison is expressed as a preference.

With preference function *P*

$$P: A \times A \to [0, 1] \tag{3}$$

the intensity of preference for alternative *a* in relation to alternative *b* is expressed, with the following interpretation:

*P* (*a, b*) = 0 marks indifference between *a* and b, that is, there is no preference of *a* over *b*,

*P* (*a, b*) *≈* 0 marks weak preference of *a* over b,

*P* (*a, b*) *≈* 1 marks strong preference of *a* over b,

*P* (*a, b*) = 1 marks strict preference of *a* over *b*.

Preference function that is added to a given criterion is the difference function of criteria value of alternatives, and it can be written as

$$P(a, b \mid) = P(k(a \mid) - k(b \mid)) = P(d \mid) \tag{4}$$

*P*(*d*) is a non-decreasing function that assumes value zero for negative difference values *d* = *k* (*a*) − *k* (*b*) , *if the functions should be maximized, that is, P* (*a, b*) = *P* (*k* (*a*) − *k* (*b*)), that is, *d* = −(*k* (*a*) − *k* (*b*)) if the criterion is minimized (**Table 4**).

<sup>5</sup> http://www.businessdictionary.com/definition/stakeholder-value-approach.html.

Weight *ti*

is the measure of relative importance of the criterion *ki*

*P* (*a, b*) *≈* 0 marks weak preference of *a* over *b* for all criteria, *P* (*a, b*) *≈* 1 marks strong preference of *a* over *b* for all criteria.

Input and output flows can be defined for every node (shown in **Figure 2**.)

Multi-criteria preference index IP is defined as the medium of preference functions *Pi*

∑ *i*=1 *k t i Pi* (*a*, *<sup>b</sup>* ) \_\_\_\_\_\_\_\_\_ ∑ *i*=1 *k t i*

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

IP (*a,b*) represents intensity, that is, the strength of decision maker's preference for activity *a* over activity *b*, when all criteria are compared at the same time. It varies between values 0 and 1.

This can also be shown in a graph. Between two nodes (two activities) *a* and *b* there are two arches with values IP(*a, b*) and IP(*b, a*). This relation is shown in **Figure 1**. There is no direct

> (*a* ) = ∑ *x*∈*k*

> (*a* ) = ∑ *x*∈*k*

*IP*(*a*, *x* )

*IP*(*a*, *x* )

value for the decision maker, all weights are equal.

*IP*(*a*, *b* ) =

connection between IP(*a, b*) and IP(*b, a*).

*T*<sup>+</sup>

*T*<sup>−</sup>

**Figure 1.** IP relation. Taken from Ref. [5].

(a) Output flow is the sum of values of output flows:

(b) Input flow is the sum of values of input flows (**Figure 3**):

Output and input flow:

. If all criteria have the same

http://dx.doi.org/10.5772/67734

:

243

**Table 4.** Types of functions in the application of the PROMETHEE1 .

#### **2.2. Multi-criteria preference index**

Let us assume that the decision maker sets preference function *Pi* and weight *t i* ; for every criterion *ki* (*i* = 1, …,*n*) of the problem (2.2).

Weight *ti* is the measure of relative importance of the criterion *ki* . If all criteria have the same value for the decision maker, all weights are equal.

Multi-criteria preference index IP is defined as the medium of preference functions *Pi* :

$$IP(a, b \mid ) = \frac{\sum\_{i=1}^{k} t\_i P\_i(a, b \mid )}{\sum\_{i=1}^{k} t\_i}$$

IP (*a,b*) represents intensity, that is, the strength of decision maker's preference for activity *a* over activity *b*, when all criteria are compared at the same time. It varies between values 0 and 1.

*P* (*a, b*) *≈* 0 marks weak preference of *a* over *b* for all criteria,

*P* (*a, b*) *≈* 1 marks strong preference of *a* over *b* for all criteria.

This can also be shown in a graph. Between two nodes (two activities) *a* and *b* there are two arches with values IP(*a, b*) and IP(*b, a*). This relation is shown in **Figure 1**. There is no direct connection between IP(*a, b*) and IP(*b, a*).

Output and input flow:

Input and output flows can be defined for every node (shown in **Figure 2**.)

(a) Output flow is the sum of values of output flows:

$$T^\*(a \mid) = \sum\_{\chi \neq k} I P(a, \chi \mid)$$

(b) Input flow is the sum of values of input flows (**Figure 3**):

$$T^{-}(a \ ) = \sum\_{x \in k} \text{IP}(a, \infty)$$

**Figure 1.** IP relation. Taken from Ref. [5].

**2.2. Multi-criteria preference index**

(*i* = 1, …,*n*) of the problem (2.2).

**Table 4.** Types of functions in the application of the PROMETHEE1

Type 6. Gauss' criterion *<sup>P</sup>*(*<sup>d</sup>* ) <sup>=</sup> <sup>1</sup> <sup>−</sup> exp{<sup>−</sup> *<sup>d</sup>*

terion *ki*

1

Taken from Ref. [5].

Let us assume that the decision maker sets preference function *Pi*

**Criterion Definition Graph**

*P*(*d* ) = {

*P*(*d* ) = {

*P*(*d* ) = {

*P*(*d* ) = ⎧ ⎪ ⎨ ⎪ ⎩

*P*(*d* ) = ⎧ ⎪ ⎨ ⎪ ⎩ \_*d <sup>m</sup>*, *<sup>d</sup>* <sup>&</sup>lt; *<sup>m</sup>* 1, *d* ≥ *m*

0, *d* ≤ *m* \_1

1, *d* > *n*

0, *d* ≤ *m* \_*<sup>d</sup>* <sup>−</sup> *<sup>m</sup> <sup>n</sup>* <sup>−</sup> *<sup>m</sup>*, *<sup>m</sup>* <sup>≤</sup> *<sup>d</sup>* <sup>≤</sup> *<sup>n</sup>*

1, *d* > *n*

.

2 \_ <sup>2</sup> *<sup>σ</sup>*<sup>2</sup>}

<sup>2</sup>, *<sup>m</sup>* <sup>&</sup>lt; *<sup>d</sup>* <sup>≤</sup> *<sup>n</sup>*

0, *d* = 0 1, *d* ≠ 0

0, *d* < *m* 1, *d* ≥ *m*

Type 1. Common criterion

242 Computer Simulation

Type 2. Quasi criterion

Type 3. Criterion with a growing linear preference

Type 4. Linear criterion with an indifference area

Type 5. Criterion with preference levels

and weight *t*

*i*

; for every cri-

• stakeholders

• weight coefficients

• suggested models.

• Potential investors (PI) • Local community (LZ).

• State (DR)

**3.1. Criteria**

• preference functions (for every criterion)

*Stakeholders* considered in ranking are as follows:

Criteria used in this study are presented in **Table 5**. These 10 criteria can be divided into two categories:

*Weight coefficients* are calculated and given in **Table 4**.

K1 Maximal usage of available potentials K2 Price according to planned installed capacity

K5 Supply safeness, expected work hours

K7 Contribution to local development and welfare

K6 Possibility of combined production of electric and thermal energy

K8 Social acceptability and sustainability of other influences on the environment

K3 Incentive purchase price K4 Technology development

K9 Period of investment return

K10 Installed power

**Table 5.** Criteria for ranking scenarios.

(2) Description criteria (K4, K6, K7 and K8).

PRMOTHEE needs criteria to be defined, according to whom the ranking will be done.

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

http://dx.doi.org/10.5772/67734

245

(1) Empirical criteria, based on the data taken from NAPOIE (K1, K2, K3, K5, K9 and K10).

Since each of the stakeholders treats each of those 10 criteria in a different manner, it is essential to define weight coefficients so that every criterion has a weight definition in relation to the stakeholder. For each of the stakeholders, the criteria were sorted into three categories:

• criteria

**Figure 2.** Output flow. Taken from Ref. [5].

**Figure 3.** Input flow. Ibid 11, Ref. [5].

## **3. Absolute ranking: PROMETHEE II**

If the decision maker wants an absolute ranking, the clear flow is considered:

Absolute ranking (PII, III) is defined in the following manner:

a PII b (a prefers b) if *T*(a) > T(b).

a III b (a is indifferent to b) if T(a) = T(b).

Elements of scientific research6 are all the elements that have to be defined so that the aforementioned mathematical model could be applied. Those comprehend:

<sup>6</sup> This research paper gave initial idea for this chapter as well as for stakeholders and used criteria [6, 7].


*Stakeholders* considered in ranking are as follows:


#### **3.1. Criteria**

**3. Absolute ranking: PROMETHEE II**

a III b (a is indifferent to b) if T(a) = T(b).

a PII b (a prefers b) if *T*(a) > T(b).

Elements of scientific research6

**Figure 3.** Input flow. Ibid 11, Ref. [5].

**Figure 2.** Output flow. Taken from Ref. [5].

244 Computer Simulation

6

If the decision maker wants an absolute ranking, the clear flow is considered:

are all the elements that have to be defined so that the afore-

Absolute ranking (PII, III) is defined in the following manner:

mentioned mathematical model could be applied. Those comprehend:

This research paper gave initial idea for this chapter as well as for stakeholders and used criteria [6, 7].

PRMOTHEE needs criteria to be defined, according to whom the ranking will be done. Criteria used in this study are presented in **Table 5**.

These 10 criteria can be divided into two categories:


*Weight coefficients* are calculated and given in **Table 4**.

Since each of the stakeholders treats each of those 10 criteria in a different manner, it is essential to define weight coefficients so that every criterion has a weight definition in relation to the stakeholder. For each of the stakeholders, the criteria were sorted into three categories:


**Table 5.** Criteria for ranking scenarios.


An assessment of weight coefficients was made on that basis, with values for K attributed on the scale of 1–10, starting from the categorization of the criteria. A representation of weight coefficients is given in **Table 6**.

**3.2. Suggested models**

2020 according to NAPOIE.

calculated and presetned in **Table 9**.

**Hydro potential**

**Total** 3360

**Table 7.** Scenarios A1–A6.

The following models (scenarios) were defined (**Table 6**):

able energy sources would be produced in mini hydros.

energy sources would be produced from biomass.

able energy sources would be produced by the Sun.

energy sources would be produced from geothermal potentials.

will be used in further researches as a continuation of this chapter.

energy sources would be produced by the wind.

• The first model (A1) represents allocation A1. This allocation fits the goals planned until

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

http://dx.doi.org/10.5772/67734

247

• The second model (A2) represents allocation A2, in which the needed energy from renew-

• The third model (A3) represents allocation A3, in which the needed energy from renewable

• The fourth model (A4) represents allocation A4, in which the needed energy from renew-

• The fifth model (A5) represents allocation A5, in which the needed energy from renewable

• The sixth model (A6) represents allocation A6, in which the needed energy from renewable

*N.B.:* It is VERY important to point out here that, according to available potentials, as shown in **Table 7** (data taken from the document 'Politika Republike Srbije u oblasti OIE'), each of the renewable energy sources listed (mini hydros, biomass, solar, wind and geothermal energy) can deliver 2252 GWh of energy independently (**Table 8** presents coneversion of available resources presented in **Table 7** from Mtoe to GWh), which represents the remainder from the total of 3360 GWh, diminished by the amount delivered by hydro potentials >10 MW. The first model A1 of this chapter was given illustratively as the goal which was set to be reached and

Scenaria are treated according to the defined criteria. Values of criteria for each scenaria are

**A1 A2 A3 A4 A5 A6 GWh GWh GWh GWh GWh GWh**

>10 MW **1108 1108 1108 1108 1108 1108** <10 MW 592 2252 0 0 0 0 **Biomass** 640 0 2252 0 0 0 **Solar** 13 0 0 2252 0 0 **Wind** 1000 0 0 0 2252 0 **Geothermal** 7 0 0 0 0 2252

*Preference functions*. A preference function is attributed to every defined criterion. Common functions according the PROMETHEE are presented in **Table 4**. For this chapter, the following allocation was adopted:



**Table 6.** Calculation of weight coefficients.

#### **3.2. Suggested models**

• Very important,

• Of little importance.

coefficients is given in **Table 6**.

ing allocation was adopted:

sible or impossible.

**State**

**Investors**

**Local community**

**Table 6.** Calculation of weight coefficients.

ference is taken as decision threshold (*m* = *d*max)

An assessment of weight coefficients was made on that basis, with values for K attributed on the scale of 1–10, starting from the categorization of the criteria. A representation of weight

*Preference functions*. A preference function is attributed to every defined criterion. Common functions according the PROMETHEE are presented in **Table 4**. For this chapter, the follow-

• Type 1. A common function is attributed to K6. Type 1 function is used when there are only two expected results, and it provides an obvious preference. Because of that it is attributed to criterion K6, since the combined production of electric and thermal energy is either pos-

• Type 3. A growing linear preference function is attributed to K2, K3, K5, K9 and K10. Type 3 function is used when the difference can be a constant value. The maximum value of dif-

• Type 4. A function with preference levels is attributed to K1, K4, K7 and K8. Type 4 function is used for discrete value differences and their outputs are discrete preferences 0, ½, 1 (*m* and *n* are decision thresholds). For criterion K1, assumed decision thresholds are *m* =

10% *d*max, and *n* = 30% *d*max, while for criteria K4, K7, K8 *m* = 1 and *n* = 2.

**Weight coefficient** *ti* **∑***ti*

k1; k5; k10 (8 + 9 + 10)/3 = 9 0.1636 Very important: 16.36%

k4; k9 (1 + 2)/2 = 1.5 0.02727 Of little importance: 2.72%

k2; k3; k4; k9 (7 + 8 + 9 + 10)/4 = 8.5 0.154545 Very important: 15.45%

k1; k7; k8 (1 + 2 + 3)/3 = 2 0.03636 Of little importance: 3.63%

k6; k7; k8 (8 + 9 + 10)/3 = 9 0.1636 Very important: 16.36% k1; k5 (6 + 7)/2 = 6.5 0.11818 Important: 11.818%

k2; k3; k4; K9; K10 (1 + 2 + 3 + 4 + 5)/5 = 3 0.0545 Of little importance: 5.45%

k5; k6; k10 (4 + 5 + 6)/3 = 5 0.0909 Important: 9.09%

k2; k3; k6; k7; k8 (3 + 4 + 5 + 6 + 7)/5 = 5 0.0909 Important: 9.09%

• Important,

246 Computer Simulation

The following models (scenarios) were defined (**Table 6**):


*N.B.:* It is VERY important to point out here that, according to available potentials, as shown in **Table 7** (data taken from the document 'Politika Republike Srbije u oblasti OIE'), each of the renewable energy sources listed (mini hydros, biomass, solar, wind and geothermal energy) can deliver 2252 GWh of energy independently (**Table 8** presents coneversion of available resources presented in **Table 7** from Mtoe to GWh), which represents the remainder from the total of 3360 GWh, diminished by the amount delivered by hydro potentials >10 MW. The first model A1 of this chapter was given illustratively as the goal which was set to be reached and will be used in further researches as a continuation of this chapter.

Scenaria are treated according to the defined criteria. Values of criteria for each scenaria are calculated and presetned in **Table 9**.


**Table 7.** Scenarios A1–A6.


**Criterion K8: social acceptability and sustainability of other influences on the environment** Most inhabitants are against any installations, regardless of their surroundings (no) 1 Inhabitants' opinion is split (split) 2 Most inhabitants accept installations, since they are far from inhabited areas and have no visible

Most inhabitants accept installations, since they are far from inhabited areas, regardless of

Mathematical model representation for the state as a stakeholder

Most inhabitants are pro installations (OK) 5

**State Min Min Min Max Max Max Max Max Min Max**

A3 Biomass 0.0861 1,591,178,750 10.74 4 6400 14 4 6.6 352 A4 Solar 0.3227 4,330,769,231 18.45 3 1300 0 1 3 10.4 1732 A5 Wind 0.9682 1,595,542,000 9.20 4 2000 0 1 3 7.7 1126 A6 Geothermal 0.9682 1,323,854,286 8.30 4 7000 1 2 5 7.1 322

**K1% K2 € K3 K4 K5 K6 K7 K8 K9 K10**

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

0.2420 1,998,203,175 9.89 5 3150 0 2 4 9.0 715

**Differentiation** *d:* **difference between scenario** *a***2 and other suggested scenarios 0.0000 0 0.00 0 0 0 0 0 0.0 0**

Biomass −0.1559 −407,024,425 0.85 1 −3250 −1 −2 0 −2.4 363 Solar 0.0807 233,266,056 8.56 2 1850 0 1 1 1.4 −1017 Wind 0.7262 −402,661,175 −0.69 1 1150 01 1 −1.3 −411 Geothermal 0.7262 −674,348,889 −1.60 1 −3850 −1 0 −1 −1.9 393

**Preference function** *P:* **scenario** *a***2 versus other suggested scenarios**

*Ti* 0.1636 0.0909 0.0909 0.02727 0.1636 0.0909 0.0909 0.0909 0.02727 0.1636

Biomass 0.0000 0 0.00 0 0 0 0 0 0.0 0 Solar 0.2366 2,739,590,481 7.71 1 5100 1 3 1 3.8 −1380

**Differentiation** *d:* **difference between scenario** *a***3 and other suggested scenarios**

**0.1559 407,024,425 −0.85 −1 3250 1 2 0 2.4 −363**

*a*3 Biomass 0 0 0.099299 0.5 0 0 0 0 0 0.923 *a*4 Solar 0 11 1 10 0.5 0.5 1 0 *a*5 Wind 1 00 0.5 0.62 0 0.5 0.5 0 0 *a*6 Geothermal 1 0 0 0.5 00 00 0 1

**0 0 0 0 00 00 0 0**

3

http://dx.doi.org/10.5772/67734

249

4

damaging effects (vis-res)

A2 Hydro potential <10 MW

*d***(***a***2,***ai***) Hydro potential <10 MW**

*P***(***a***2,***ai***)Hydro potential <10 MW**

*d***(***a***3,***ai***)Hydro potential <10 MW**

whether there is a visual contact (res)

**Table 8.** Available potentials of renewable energy sources.


**Table 9.** Scenarios according to K criteria values.

## **4. Mathematical model**



#### Mathematical model representation for the state as a stakeholder


#### *d***(***a***2,***ai***) Hydro potential Differentiation** *d:* **difference between scenario** *a***2 and other suggested scenarios**


#### *P***(***a***2,***ai***)Hydro potential Preference function** *P:* **scenario** *a***2 versus other suggested scenarios**

**4. Mathematical model**

**Table 9.** Scenarios according to K criteria values.

companies in energy production sector)

industrial regions on wider areas)

improvement)

248 Computer Simulation

market (com\_EU)

(moderate)

**Criterion K4: Technology development**

A1 PLAN 43.00 1,356,627,968 9.87 4 3564 1 34 6.1 799 A2 Hydro potential <10 MW 24.20 1,998,203,175 9.89 5 3150 0 2 4 9.0 715 A3 Biomass 8.61 1,591,178,750 10.74 4 6400 1 44 6.6 352 A4 Solar 32.27 4,330,769,231 18.45 3 1300 0 1 3 10.4 1732 A5 Wind 96.82 1,595,542,000 9.20 4 2000 0 1 3 7.7 1126 A6 Geothermal 96.82 1,323,854,286 8.30 4 7000 1 2 5 7.1 322

**K1 (%) K2 (€) K3 K4 K5 K6 K7 K8 K9 K10**

**Criterion K7: Contribution to local development**

3

5

3

4

5

Technologies in laboratory and research phases (laboratory) 1 Technologies in pilot programs (pilot) 2

Commercially ready technologies with a reliable place in the overall local market (com\_loc) 4

Without any influence on local economy (none) 1 Weak influence on local economy(weak) 2

Commercially ready technologies with a reliable place in the supranational and European

Moderate influence on local economy (only a small number of permanent workplaces)

Moderate to large influence on local economy (opening new workplaces and chains of

Very large influence on local economy (strong incentive to local growth, creation of small

Technologies demanding further improvements to enhance their efficiency (further

**Type of renewable energy sources Mtoe GWh** Hydro 0.8 9304 Biomass 2.25 26,167 Solar 0.6 6978 Wind 0.2 2326 Geothermal energy 0.2 2326

**Table 8.** Available potentials of renewable energy sources.


#### *d***(***a***3,***ai***)Hydro potential Differentiation** *d:* **difference between scenario** *a***3 and other suggested scenarios**



*P***(***a***3,***ai***) Hydro potential** 

*d***(***a***6,***ai***) Hydro potential <10 MW**

*P***(***a***6,***ai***) Hydro potential <10 MW**

**N.B.:** IP = (*ai,as*), *i,s* = 2,3,4,5,6; IP = ∑*tjPj*(*ai,as*).

**<10 MW**

0.02727 0.1636 0.0909 0.0909 0.0909 0.02727 0.1636

**Preference function** *P* **– scenario** *a***5 versus other suggested scenarios**

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

http://dx.doi.org/10.5772/67734

251

0 0.147 0.075 0 0 0 0 0 0.481 0.51 *a*3

**Differentiation** *d* **– difference between scenario** *a***6 and other suggested scenarios**

Biomass −0.8821 267,324,464 2.45 0 600 0 −2 1 −0.5 −30 Solar −0.6455 3,006,914,945 10.16 1 5700 1 1 2 3.3 −1410 Wind 0.0000 271,687,714 0.90 0 5000 11 2 0.6 −804 Geothermal 0.0000 0 0.00 0 0 0 0 0 0.0 0

**Preference function** *P* **– scenario** *a***6 versus other suggested scenarios**

*ti* 0.1636 0.0909 0.0909 0.02727 0.1636 0.0909 0.0909 0.0909 0.02727 0.1636

*a*3 Biomass 0 0.089 0.241 0 0.105 0 0 0.5 0 0 *a*4 Solar 011 0.5 1 1 0.5 11 0 *a*5 Wind 0 0.09 0.088 0 0.877 1 0.5 1 0.182 0 *a*6 Geothermal 0 0 0 0 0 0 0 0 0 0

**IP(***a***2,***a***3) IP(***a***2,***a***4) IP(***a***2,***a***5) IP(***a***2,***a***6)** 0.1736 0.49084 0.369567 0.340835 IP(*a*3,*a*2) IP(*a*3,*a*4) IP(*a*3,*a*5) IP(*a*3,*a*6) 0.398447 0.695355 0.5399 0.421672 IP(*a*4,*a*2) IP(*a*4,*a*3) IP(*a*4,*a*5) IP(*a*4,*a*6) 0.117956 0.160001 0.233948 0.3272 IP(*a*5, *a*2) IP(*a*5,*a*3) IP(*a*5,*a*4) IP(*a*5,*a*6) 0.116733 0.172473 0.386305 0.1636 IP(*a*6,*a*2) IP(*a*6,*a*3) IP(*a*6,*a*4) IP(*a*6,*a*5) 0.29712 0.92625 0.613555 0.391871

−0.7262 674,348,889 1.60 −1 3850 1 0 1 1.9 −393

0 0.224 0.157 0 0.675 1 0 0.5 0.576 0

Biomass 0 0 0.166 0 0 0 0 00 0.962 *a*4 Solar 0 11 0.5 1 0 0 0 10 *a*5 Wind 0 0 00 0 0 0 0 00 *a*6 Geothermal 0 0 0 00 000 01 *ti* 0.1636 0.0909 0.0909

#### *P***(***a***3,***ai***) Hydro potential Preference function** *P***: scenario** *a***3 versus other suggested scenarios**






*d***(***a***6,***ai***) Hydro potential <10 MW Differentiation** *d* **– difference between scenario** *a***6 and other suggested scenarios**





**N.B.:** IP = (*ai,as*), *i,s* = 2,3,4,5,6; IP = ∑*tjPj*(*ai,as*).

*d***(***a***3,***ai***)Hydro potential <10 MW**

250 Computer Simulation

*P***(***a***3,***ai***) Hydro potential <10 MW**

*d***(***a***4,***ai***) Hydro** 

**potential<10 MW**

*P***(***a***4,***ai***) Hydro potential <10 MW**

*d***(***a***5,***ai***) Hydro potential <10 MW**

**Differentiation** *d:* **difference between scenario** *a***3 and other suggested scenarios**

Wind 0.8821 4,363,250 −1.54 0 4400 1 3 1 1.1 −774 Geothermal 0.8821 −267,324,464 −2.45 0 −600 0 2 −1 0.5 30

**Preference function** *P***: scenario** *a***3 versus other suggested scenarios**

*ti* 0.1636 0.0909 0.0909 0.02727 0.1636 0.0909 0.0909 0.0909 0.02727 0.1636

Biomass −0.2366 −2,739,590,481 −7.71 −1 −5100 −1 −3 −1 −3.8 1380

Solar 0.0000 0 0.00 0 0 0 0 0 0.0 0 Wind 0.6455 −2,735,227,231 −9.25 −1 −700 00 0 −2.7 606 Geothermal 0.6455 −3,006,914,945 −10.16 −1 −5700 −1 −1 −2 −3.3 1410

**Preference function** *P***: scenario** *a***4 versus other suggested scenarios**

*a*3 Biomass 0 0 0 0 0 0 0 0 0 0.978

*a*4 Solar 0 00 0 0 0 00 0 0

*a*6 Geothermal 1t 0 0 0 0 0 0 0 0 1

*a*5 Wind 1 0 0 0 0 0 0 0 0 0.43

*ti* 0.1636 0.0909 0.0909 0.02727 0.1636 0.0909 0.0909 0.0909 0.02727 0.1636

Biomass −0.8821 −4,363,250 1.54 0 −4400 −1 −3 −1 −1.1 774 Solar −0.6455 2,735,227,231 9.25 1 700 0 0 0 2.7 −606

Wind 0.0000 0 0.00 0 0 0 0 0 0.0 0 Geothermal 0.0000 −271,687,714 −0.90 0 −5000 −1 −1 −2 −0.6 804

**Differentiation** *d***: difference between scenario** *a***5 and other suggested scenarios −0.7262 402,661,175 0.69 −1 −1150 0 −1 −1 1.3 411**

*a*3 Biomass 0 0 0 0 0 0 0 0 0 0 *a*4 Solar 0.5 1 1 0.5 1 11 0.5 1 0 *a*5 Wind 1 0.0016 0 0 0.862 1 1 0.5 0.289 0 *a*6 Geothermal 1 0 0 0 0 0 1 0 0.131 1

**0.1559 407,024,425 −0.85 −1 3250 1 2 0 2.4 −363**

**0.5 0.148 0 0 0.637 1 1 0 0.63 0**

**Differentiation** *d***: difference between scenario** *a***4 and other suggested scenarios −0.0807 −2,332,566,056 −8.56 −2 −1850 0 −1 −1 −1.4 1017**

**0 0 0 0 0 0 0 0 0 0.721**


*a***2** *a***3** *a***4** *a***5** *a***6 T+ T** *a*6 0.294074 0.041479 0.622703 0.267196 0 0.306363 0.000179

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

http://dx.doi.org/10.5772/67734

253

T- 0.2526075 0.124051 0.582936 0.289638 0.306185

**Figure 5.** Chart representation of ranking results for the investors as stakeholders.

IP(*a*2,*a*3) IP(*a*2,*a*4) IP(*a*2,*a*5) IP(*a*2,*a*6) 0.082911 0.393418 0.316357 0.19993 IP(*a*3,*a*2) IP(*a*3,*a*4) IP(*a*3,*a*5) IP(*a*3,*a*6) 0.436219 0.670658 0.553205 0.34342 IP(*a*4,*a*2) IP(*a*4,*a*3) IP(*a*4,*a*5) IP(*a*4,*a*6) 0.039295 0.053301 0.141615 0.17268 IP(*a*5,*a*2) IP(*a*5,*a*3) IP(*a*5,*a*4) IP(*a*5,*a*6) 0.066109 0.061476 0.202568 0.0545 IP(*a*6,*a*2) IP(*a*6,*a*3) IP(*a*6,*a*4) IP(*a*6,*a*5) 0.305534 0.10103 0.611568 0.438984

*a***2** *a***3** *a***4** *a***5** *a***6 T+ T** *a*2 0 0.082911 0.393418 0.316357 0.19993 0.248154 0.036365 *a*3 0.436219 0 0.670658 0.553205 0.34342 0.500876 0.426196 *a*4 0.039295 0.053301 0 0.141615 0.17268 0.101723 −0.36783 *a*5 0.066109 0.061476 0.202568 0 0.0545 0.096163 −0.26638 *a*6 0.305534 0.10103 0.611568 0.438984 0 0.364279 0.171647

T− 0.211789 0.0746795 0.469553 0.36254 0.192633

**Figure 6.** Chart representation of ranking results for the local community as a stakeholder.

The results for community are shown in **Figure 6**.

The reuslts for investors are shown in **Figure 5**.

*For the local community as a stakeholder:*

**Determination of preference index**

IP = (*ai, as*), *i,s* = 2,3,4,5,6; IP= ∑*tjPj*(*ai, as*).

The results for State are shown in **Figure 4**.

**Figure 4.** Chart representation of ranking results for the state as a stakeholder.

The same approach could be usd for detailed calculation for the investors and local community as stakeholders.


*For the investors as stakeholders:*

N.B.: IP = (*ai, as*), *i, s* = 2,3,4,5,6; IP= ∑*tjPj*(*ai,as*).



The reuslts for investors are shown in **Figure 5**.

*a***2** *a***3** *a***4** *a***5** *a***6 T+ T** *a*2 0 0.1736 0.49084 0.369567 0.340835 0.343711 0.111147 *a*3 0.398447 0 0.695355 0.5399 0.421672 0.513844 0.364169 *a*4 0.117956 0.160001 0 0.233948 0.3272 0.209776 −0.33674 *a*5 0.116733 0.172473 0.386305 0 0.1636 0.209778 −0.17404 *a*6 0.29712 0.092625 0.613555 0.391871 0 0.348793 0.035466

The same approach could be usd for detailed calculation for the investors and local commu-

*a***2** *a***3** *a***4** *a***5** *a***6 T+ T**

*a*2 0 0.1611 0.590895 0.272998 0.359078 0.346018 0.09341 *a*3 0.377197 0 0.640883 0.402209 0.33841 0.439675 0.315624 *a*4 0.195746 0.206178 0 0.21615 0.281805 0.22497 −0.35797 *a*5 0.143413 0.087446 0.477263 0 0.245445 0.238392 −0.05125

IP(*a*2,*a*3) IP(*a*2,*a*4) IP(*a*2,*a*5) IP(*a*2,*a*6) 0.1611 0.590895 0.272998 0.359078 IP(*a*3,*a*2) IP(*a*3,*a*4) IP(*a*3,*a*5) IP(*a*3,a6) 0.377197 0.640883 0.402209 0.33841 IP(*a*4,*a*2) IP(*a*4,*a*3) IP(*a*4,*a*5) IP(*a*4,*a*6) 0.195746 0.206178 0.21615 0.281805 IP(*a*5,*a*2) IP(*a*5,*a*3) IP(*a*5,*a*4) IP(*a*5,*a*6) 0.143413 0.087446 0.477263 0.245445 IP(*a*6,*a*2) IP(*a*6,*a*3) IP(*a*6,*a*4) IP(*a*6,*a*5) 0.294074 0.041479 0.622703 0.267196

T− 0.232564 0.149675 0.546514 0.383822 0.313327

**Figure 4.** Chart representation of ranking results for the state as a stakeholder.

The results for State are shown in **Figure 4**.

nity as stakeholders.

252 Computer Simulation

*For the investors as stakeholders:* **Determination of preference index**

N.B.: IP = (*ai, as*), *i, s* = 2,3,4,5,6; IP= ∑*tjPj*(*ai,as*).

**Figure 5.** Chart representation of ranking results for the investors as stakeholders.


*For the local community as a stakeholder:*

The results for community are shown in **Figure 6**.

**Figure 6.** Chart representation of ranking results for the local community as a stakeholder.

## **5. Chart representation of results**

After applying the PROMETHEE, as a tool for stakeholder value approach, and after the ranking, we can reach following conclusions on the basis of results obtained:

**References**

r-edward-freeman/

ljivih- izvora-energije

approach.html

Nauka, Univerziteta u Beogradu; 2010.

Policy. 2009;**37**(2009):1587–1600.

[1] Edward Freeman, R. Univesity of Virginia. Faculty & Research. The Faculty Directory [Internet]. Available from: http://www.darden.virginia.edu/faculty-research/directory/

Multi-Criteria Decision-Making in the Implementation of Renewable Energy Sources

http://dx.doi.org/10.5772/67734

255

[2] Radna Grupa.Ministarstvo energetike rudarstva i zastite zivotne sredine. Nacionalni akcioni plan za primenu obnovljivih izvora energije [Internet]. June 2013. Available from:

[3] Narodna Skupština Republike Srbije zakon o ratifikaciji. Zakon o ratifikaciji ugovora o osnivanju energetske zajednice i Republike Albanije, Rapublike Bugarske, Bosne i Hercegovine, Republike Hrvatske, Bivše Jugoslovenske Republike Makednoije, Republike Crne Gore, Rumunije, Republike Srbije, i privremene misije Ujedinjenih Nacija na Kosovu u skladu sa rezolucijom 1244 Saveta bezbednosti Ujedinjenih Nacija [Internet]. July 2006. Available

[4] Zorana Mihajlovic. Mistarstvo energetike rudarstva i zastite zivotne sredine. Politika Republike Srbije u oblasti obnovljivih izvora energije [Internet]. December 2012. Available from: http://www.serbio.rs/item/102-politika-republike-srbije-u-oblasti-obnov-

[5] Milutin Čupić, Milija Suknović. Odlučivanje. 6th ed. Beograd: Fakultet organizacionih

[6] Available from: http://www.businessdictionary.com/definition/stakeholder-value-

[7] T. Tsoutsos T, Drandaki M, Frantzeskaki N, Iosifidis E, Kiosses I. Sustainable energy planning by using multi-criteria analysis application in the island of Crete. Energy

http://www.mre.gov.rs/dokumenta-efikasnost-izvori.php

from: http://www.parlament.gov.rs/narodna-skupština. 1.html

The results obtained and shown in the charts indicate, in fact, that, according to defined criteria and weight coefficients, the sequence of types of renewable energy sources is absolutely identical regardless of the stakeholder. The sequence of priorities in the application of renewable energy sources for the production of electricity goes as follows:


Further activities of all stakeholders should be given to mini hydros and biomass, since they have the best relation toward the aforementioned criteria.

According to presented model, potentials of all the mentioned types of renewable energy sources are capable for achieving its goals, with the limitation that wind and geothermal energy would have, according to such a premise, a 96.82% usage, which is not a convenient circumstance, while biomass would have an 8.61% usage and mini hydros 24.20%.

The general conclusion is that the state as a stakeholder should focus its activities regarding the production of electricity from renewable energy sources on biomass and mini hydros, since, according to listed hypotheses, defined criteria and the application of the mathematical model, they proved to be the best solution. The same goes for investors and local community as stakeholders.

Methodology use in this chapter is taken into account the criteria and stakeholders which where possible to use according to the official available data. The final number of stakolders and criteria are endless and just make calculation model more comprehensive.

## **Author details**

Dejan Jovanovic1 \* and Ivan Pribicevic2

\*Address all correspondence to: deki.jovanovic@gmail.com

1 JP Zavod za udzbenike, Belgrade, Serbia

2 TMS CEE d.o.o Beograd, Belgrade, Serbia

## **References**

**5. Chart representation of results**

(1) Biomass

254 Computer Simulation

(2) Mini hydros (3) Geothermal (4) Wind energy (5) Solar energy

as stakeholders.

**Author details**

Dejan Jovanovic1

After applying the PROMETHEE, as a tool for stakeholder value approach, and after the rank-

The results obtained and shown in the charts indicate, in fact, that, according to defined criteria and weight coefficients, the sequence of types of renewable energy sources is absolutely identical regardless of the stakeholder. The sequence of priorities in the application of renew-

Further activities of all stakeholders should be given to mini hydros and biomass, since they

According to presented model, potentials of all the mentioned types of renewable energy sources are capable for achieving its goals, with the limitation that wind and geothermal energy would have, according to such a premise, a 96.82% usage, which is not a convenient

The general conclusion is that the state as a stakeholder should focus its activities regarding the production of electricity from renewable energy sources on biomass and mini hydros, since, according to listed hypotheses, defined criteria and the application of the mathematical model, they proved to be the best solution. The same goes for investors and local community

Methodology use in this chapter is taken into account the criteria and stakeholders which where possible to use according to the official available data. The final number of stakolders

circumstance, while biomass would have an 8.61% usage and mini hydros 24.20%.

and criteria are endless and just make calculation model more comprehensive.

ing, we can reach following conclusions on the basis of results obtained:

able energy sources for the production of electricity goes as follows:

have the best relation toward the aforementioned criteria.

\* and Ivan Pribicevic2

1 JP Zavod za udzbenike, Belgrade, Serbia 2 TMS CEE d.o.o Beograd, Belgrade, Serbia

\*Address all correspondence to: deki.jovanovic@gmail.com


## *Edited by Dragan Cvetković*

The first chapter provides an overview of the development of a novel agent-based simulation model of socio-environmental innovation diffusion. The second chapter shows the study about rendering of colours with three rendering engines. The third and fourth chapters are devoted to modelling clothes at different levels. The fifth chapter describes the modelling of computer simulation in the optimization of bioprocess technology. Chapters 6 and 7 formulate a physical model of deformation of steel and idea of constructing a scientific workshop focused on high-temperature processes. Chapter 8 formulates surrogate models. Chapter 9 shows computer simulation of high-frequency electromagnetic fields. Chapter 10 proposes the modelling of the task allocation problem by the use of Petri Nets. Chapter 11 presents various scenarios whose ranking is done according to defined criteria and weight coefficients.

Photo by gleitfrosch / iStock

Computer Simulation

Computer Simulation