**2. Indonesia**

Human echinococcosis in Indonesia was first reported in nonindigenous people, in a European with liver hydatidosis in 1885 [7]. Subsequent cases were also confirmed in four European and an Indian in the 1920s [8]. In the late 1930s, a 45-yearold European female was also suspected with *Echinococcus* infection after stayed for almost 17 years which much contact with dogs in Indonesia [9].

The first indigenous case was diagnosed in a female aged 23 years from the northern part of South Sulawesi in 1988. She was admitted to the hospital with pain in the right under quadrant of the abdomens and in the right part of the hypochondrium, below the right rib. Hepatomegaly or tumor in the right upper quadrant of the abdomen was not found by physical examination. However, a hypoechoic, oval image, with a clear acoustic image on the right, was identified by using ultrasonography. Serological test was not performed, and the diagnosis was confirmed by radiological image examination. The patient has a history of keeping a dog [10].

Another echinococcosis case was reported in a 90-year-old female from Java in 1988. She was hospitalized in Ujung Pandang (Makassar), South Sulawesi with complaint of severe weakness, fainting, and no appetite. A tumor with no malignant cells was identified in the lower right part of the abdomen. Cysts with hypoechoic image on the walls were detected in the pelvis, right upper liver, and upper left abdomen. Fluid aspiration and serological test were not carried out during patient treatment at the hospital. The hydatidosis was confirmed in the right diaphragm arc/liver, omentum, and adnexa by radio-image. She had not kept dog at home [11].

In 1997, subsequent indigenous case, a 26-year-old female from Palu, Central Sulawesi, was administered in the Cipto Mangunkusumo Hospital in Jakarta with chronic and progressive headache, space-occupying lesion, blurred vision since the last 6 months, and vomiting. Bilateral edema of the papilla of the eyes was found without other neurological deficits. Cysts were found in many parts of organs. A hypodense, cystic lesion in the left parietal region of the brain divided by several septa with calcifications was discovered by radio-image with tomography (CT) scan. Historically, he lived in an area that kept many dogs, including her house.

In animals, echinococcosis is reported from goat, sheep, pig, water buffalo, and monkey from Sumatra, Madura, Bali, and Sulawesi. A case with liver, lung, and spleen hydatidosis was first identified in a cattle in Buitenzorg (Bogor, West Java) during meat inspection on 1906. Hydatid cysts were subsequently found from the lungs and liver of cattle in North Sumatra in 1907 [12]. At that time, the hydatid cysts were only detected from imported Australian cattles (nonindigenous echinococcosis was also reported among 46 cattles and 3 pigs in Medan, North Sumatra, in 1916) [13]. In 1929, hydatid cysts in indigenous cattle were found for the first time from Bali. This finding was considered that dogs from Bali were probably reservoirs for the adult worm *Echinococcus granulosus.* However, the lack of study has caused this to never to be proven [14]. Further cases of echinococcosis in animals were reported back in 1940. A single cyst of *Echinococcus* was found in the lung of the goat (*Capra hircus*) and the lung of the buffalo (*Buffelus sondaicus*) in the

Netherlands Indies [15]. In 1947, a total of 45 cysts of *Echinococcus unilocularis* (*E. granulosus*) were found in the lungs, kidney, stomach wall, and liver on an old monkey *Cynopithecus niger* (later on as *Macaca tonkeana*) from Sulawesi (Celebes) that suddenly died. The first evidence of *E. granulosus* infection in domestic dogs in Indonesia was reported from the Lindu Valley, Central Sulawesi, in 1972. The dog died suddenly after 18 months with numerous necropsies due to *E. granulosus.* Subsequent dog with echinococcosis was reported on the shore of Lake Lindu at Tomado village in 1973. The positive eggs of *E. granulosus* were confirmed by direct smear and formaldehyde-ether from the dog stools. The adult worms of *E. granulosus* were also identified during the study. The worms were collected from a dead dog on the shore of Lake Lindu at Tomado [16]. Since then, hundreds of dogs have been examined to identify the infection rate of *E. granulosus* in this area. However, hydatid cysts and adult *E. granulosus* have not been discovered until now [13].


**4. Building Human Capabilities**

Research on capabilities can inform public policy. This special issue of *Education Science* focuses on research that informs social action to build human capabilities for people left behind in their local context and promotes learning across contexts.

Human capabilities emerged from comparisons of states in India and women's opportunities in Islamic countries [112–114]. Nussbaum argued that women had a right to education to a level needed to support a family, a threshold of wellbeing not met in countries that denied education to women. Nussbaum [112,113] focuses on health, education, faith, and other freedoms and responsibilities. Building on her ideas, the capabilities noted in Table 3 relate to education and career pathways, along with localized social support for actualizing them. Local schools and community organizations are at the center of this approach. Their capacity to organize interventions is enhanced by investment from the government, businesses, and foundations, including university–school partnerships that provide information and local groups providing mentoring and social support (Table 3).

**Table 3.** Human-capabilities framework for addressing educational inequalities: assessing gaps and opportunities for social support and networking.

Comparing strategies used to support marginalized students and communities in Ireland and the US [61] reveals no simple solutions to K-12 preparation and access inequality. Of course, the US is much larger but about half of the states have more population than Ireland. A comparison of US states also shows that states that coordinated student aid funding with rising costs and provided social support saw more substantial gains in college enrollment rates than comparable states [43,114].

The human-capabilities approach provides workable means for research informing policy and practice. Studies informing policymakers in Ireland, England, and the US confirm the workability of this approach in raising capabilities (preparation, access, and college success) for groups that would likely otherwise be left behind [9,69,115]. Students' engagement in encouragement and support services improves the odds of college enrollment [47]. Research informs practitioners about what does not work as intended, why this is the case, and alternative policies and practices. However, substantial changes will likely be needed to transport innovations across contexts, even when successful interventions occur.

In addition, and perhaps more importantly, action studies informing social support projects help expand and refine services and the quality of information and mentoring they provide. The global experience of the move to marketization through student loans in higher education raises differing social issues across nations. Consider the example of possible social interventions that would minimize the negative consequences of growing student debt, a problem of emerging importance across countries in the Pacific region. The economic wellbeing of students who graduate with debt, as measured by their ability to support their families, is a primary indicator in the human-capabilities framework (Table 3). Since local and national economies vary substantially, the local remedies to this challenge must differ to fit policy circumstances and cultures.


**4.2 Transient flow**

In this section the transient simulations results are presented. **Figure 7** represents the evolution of vapor volume fraction at different times for L2/d2 = 5. The inlet corner radius rc is equal to zero. The injection pressure *pin* = 1000 bar and the exit pressure *pout* = 50 bar.

The cavitation appears in the vicinity of the sharp edge for a time of the order of 0.6 μs. Then, the cavitation pocket elongates progressively through the orifice and reaches the nozzle exit at t = 3 μs. This result agrees well with the numerical

*Analysis of Geometric Parameters of the Nozzle Orifice on Cavitating Flow and Entropy… DOI: http://dx.doi.org/10.5772/intechopen.99404*

**Figure 6.** *Discharge coefficient as a function of rc/d2 ration for Re = 2 103 and Re = 6.8 103 .*

**Figure 7.** *Transient evolution of vapor volume fraction.*

simulation obtained by Dumont et al. [37] in a similar injector and by experimental visualization and measurements Ohrn et al. [33].

**Figure 8** depicts the axial profile of the mixture density near the wall at t = 1.5 μs and t = 3.8 μs. Upstream of the corner the fluid is at liquid state.

At the sharp edge, the density exhibits an important decrease due to the pressure drop and the appearance of cavitation. In this region, the mixture density contains mainly fuel vapor.

Downstream the corner, the mixture density increases and the vapor fraction decreases, owing to the collapse and the breakup of the cavitation along the orifice.

At the nozzle exit, the cavitation is present for t = 3.8 μs and therefore the mixture density is smaller than the liquid density *ρl*. Thus the spray that leaves the injector and penetrates into the combustion chamber is formed not only by fuel liquid but also by fuel vapor.

**Figure 8.** *Mixture density profile at t = 1.5 μs and t = 3.8 μs near the nozzle wall.*


Section 1 Basic Knowledge

**3**

**Chapter 1**

**Abstract**

**1. History**

**1.1 Introduction**

History and Genetics of

*Leen Abu Safieh and Saad Aldahmash*

and counseling for patients with retinoblastoma.

epigenetics, Knudson, allele, germline, somatic

*Tariq Alzahem, Waleed Alsarhani, Abdullah Albahlal,*

The history of retinoblastoma (RB) goes back to 1597 when Pieter Pawius of Amsterdam described a tumor that resembled retinoblastoma. "Fungus haematodes" was the first term used to describe retinoblastoma. Later, the American Ophthalmological Society approved the term retinoblastoma in 1926. The retinoblastoma protein is encoded by the RB1 gene located at 13q14. The functioning model of the tumor suppressor genes was first proposed by Alfred Knudson in the 1970s who precisely explained the hereditary mechanism of retinoblastoma. If both alleles of this gene are mutated, the protein is inactivated and this results in the development of retinoblastoma. One mutation can be either germline or somatic and the second one is always somatic. Differentiation between sporadic and germline retinoblastoma variants requires the identification of the RB1 germline status of the patient. This identification is important for assessing the risk of additional tumors in the same eye, the other eye, and the risk of secondary tumors. Thus, genetic testing is an important component of the management of all children diagnosed with retinoblastoma. In this chapter, we will go over the history, genetics,

**Keywords:** history, genetics, two-hit hypothesis, RB1 gene, tumor suppressor,

Different tumors that resembled retinoblastoma were described in the past, the first of which was by Pieter Pawius of Amsterdam in 1597 [1]. The Dutch anatomist described a malignancy involving the left eye, temporal bone, and cranium [1]. Pawius described the tumor as "filled with substance like brain tissue mixed with thick blood and like crushed stone." The crushed stone may have represented retinoblastoma classic calcifications. Then, in 1767, Hayes described another case of a child who had presented with the appearance "cat's eyes in the dark" [2]. This was the first description of leukocoria in the scientific literature. In the last two centuries, retinoblastoma had different names. "Fungus haematodes" was the term used in the early 1800s. This term was used to describe a fungating vascular tumor that affected different parts of the body including the globe. Before the era of ophthalmoscopy, no one thought the tumor originated from the globe except

Retinoblastoma


*Article* **Mapping of Pro-Equity Interventions Proposed by Immunisation Programs in Gavi Health Systems Strengthening Grants**

**Joelle Ducharme 1,\*, Gustavo Caetano Correa 1, Heidi W. Reynolds 1, Alyssa B. Sharkey 2, Virginia A. Fonner 3 and Mira Johri 4,5**

**Abstract:** Reaching zero-dose (ZD) children, operationally defined as children who have not received a first dose of the diphtheria, tetanus, and pertussis (DTP1) vaccine, is crucial to increase equitable immunisation coverage and access to primary health care. However, little is known about the approaches already taken by countries to improve immunisation equity. We reviewed all Health System Strengthening (HSS) proposals submitted by Gavi-supported countries from 2014 to 2021 inclusively and extracted information on interventions favouring equity. Pro-equity interventions were mapped to an analytical framework representing Gavi 5.0 programmatic guidance on reaching ZD children and missed communities. Data from keyword searches and manual screening were extracted into an Excel database. Open format responses were analysed using inductive and deductive thematic coding. Data analysis was conducted using Excel and R. Of the 56 proposals included, 51 (91%) included at least one pro-equity intervention. The most common interventions were conducting outreach sessions, tailoring the location of service delivery, and partnerships. Many proposals had "bundles" of interventions, most often involving outreach, microplanning and community-level education activities. Nearly half prioritised remote-rural areas and only 30% addressed genderrelated barriers to immunisation. The findings can help identify specific interventions on which to focus future evidence syntheses, case studies and implementation research and inform discussions on what may or may not need to change to better reach ZD children and missed communities moving forward.

**Keywords:** immunization; vaccines; zero-dose; equity; health systems


**4. Phase domain line model**

274 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

**Figure 2.** Frequency domain circuit representation of a multi-conductor line.

instance, the transformation of (21) to the time domain yields:

represented as follows:

auxiliary sources *Iaux,0* and *Iaux,L*.

currents can be represented as

with

and

with, *Ish,0* =*YCV0* , *Iaux,0* =*HIrfl,L* ,and *Irfl,L*= *IL+ YCVL* .

In the same way as it has been previously done for (19), expression (20) is conveniently

0 ,0 ,0 *sh aux* **II I** (22)

Expressions (21) and (22) constitute a traveling wave line model for the segment of length *L* depicted in figure 1. The former set of expressions represents end *L* of the line segment, while the latter set represents end *0*. A schematic representation for the whole model is provided by figure 2. Note that the coupling between the two line ends is through the

The line model defined by expressions (21) and (22) is in the frequency domain. Power system transient simulations require this model to be transformed to the time domain. For

Note that at (23), (24), (25) the lowercase variables represent the time domain images of their uppercase counterparts at (22) and that the symbol \* represents convolution. Reflected

Expressions (23)-(26) constitute a general traveling–wave based time–domain model for line end *0*. The model corresponding to the other end is obtained by interchanging sub–indexes "*0*" and "*L*" at (23)-(26). Equation (23) essentially provides the interface of the line–end *0* model to the nodal network solver that, for power system transient analysis, usually is the EMTP (Dommel, 1996). Expressions (24) and (25) require the performing of matrix–to–vector

0 ,0 ,0 *sh aux* **ii i** (23)

,0 \* *sh C L* **i yv** (24)

,0 ,0 \* *aux rfl* **i hi** (25)

, , , 2 *rfl L sh L aux L* **i ii** (26)

Since rational fitting and model solutions are carried out directly in the phase domain, the model described here is said to be a phase domain line model. Rational fitting for this model is carried out using the Vector fitting (VF) tool (Gustavsen, 2008). In the case of *YC,* the whole fitting process is done in the phase domain, whereas for *H* initial poles and time delays are first calculated in the modal domain.

#### **4.1. Rational approximation of Yc**

The following rational representation has been proposed for *YC* in (Marti, 1982) and (Morched et al., 1999):

$$\mathbf{Y}\_{\mathcal{C}} = \mathbf{G}\_0 + \sum\_{i=1}^{Ny} \frac{\mathbf{G}\_i}{\mathbf{s} - \boldsymbol{\eta}\_i} \tag{27}$$

where *Ny* is the fitting order, *qi* represents the *i–th* fitting pole, *Gi* is the corresponding matrix of residues and *G0* is a constant matrix obtained at the limit of *YC* when *s=j*. Note in (27) that common poles are used for the fitting of all the elements of *YC* obtained by fitting the matrix trace and finally the fitting of the matrices of residues *Gi* and of proportional terms *G0* is done in the phase domain. Section 7 provides an overview of the VF procedure and further information is to be found in (Gustavsen & Semlyen, 1999) and (Gustavsen, 2008).

#### **4.2. Rational approximation of H**

Rational fitting of transfer matrix *H* is substantially more involving than the one of *YC* above. To attain an accurate and compact (low order) rational representation for *H* it is essential to factor out all terms involving time delays (Marti, 1982). The major difficulty here is that its elements could involve a mix of up to *N* different delay terms due to the multimode propagation on an *N*–conductor line (Wedepohl, 1965). Separation of matrix *H* into singledelay terms is obtained from the following modal factorization (Wedepohl, 1965):

$$\mathbf{H} = \mathbf{M} \mathbf{H}\_m \mathbf{M}^{-1} \tag{28}$$

where *Hm* is a diagonal matrix of the form

$$\mathbf{H}\_m = \text{diag}\{e^{-\gamma\_1 L}, e^{-\gamma\_2 L}, \dots, e^{-\gamma\_N L}\} \tag{29}$$

and =(*YZ*) is the propagation constant of a conductor line (Wedepohl, 1965). As the triple product in (28) is performed by partitioning *M* in columns and *M–1* in rows, the following expression is derived:

$$\mathbf{H} = \sum\_{i=1}^{N} \mathbf{D}\_i e^{-\boldsymbol{\gamma}\_i L} \tag{30}$$

where *Di* is the rank–1 matrix obtained from pre–multiplying the *i–th* column of *M* by the *i– th* row of *M–1*. Matrix *Di* in fact is an idempotent (Marcano & Marti, 1997). The exponential factors at (30) can be further decomposed as follows:

$$e^{-\gamma\_i L} = e^{-\gamma\_i L} e^{-s\tau\_i}; i = 1, 2, \dots, N \tag{31}$$

where exp *i L* is a minimum phase–shift function (Bode, 1945) and *i* is the time delay associated to the velocity of the *i–th* mode. Hence:

$$\mathbf{H} = \sum\_{i=1}^{N} \mathbf{D}\_i e^{-\mathcal{F}\_i L} e^{-s\mathbf{r}\_i} \tag{32}$$

The time delays in (32) can be initially estimated by applying Bode's relation for minimum phase complex functions (Bode, 1945) to the magnitudes of *exp(-iL)* in (30). Although (32) provides the desired separation of *H* as a sum of terms, each one involving a single delay factor, the following consideration is brought in for computational efficiency (Morched et al., 1999). Modal delays often occur in groups with almost identical values. Suppose that a number *Ng* of these groups can be formed and that (32) can be modified as follows:

$$\mathbf{H} = \sum\_{k=1}^{N\_{\mathcal{S}}} \tilde{\mathbf{H}}\_k e^{-s\tau\_k} \tag{33}$$

where *Ng* is less or equal to *N*, and *k* is the representative delay for the *k–th* group. Clearly, by comparing (33) and (32):

$$\tilde{\mathbf{H}}\_k = \sum\_{i=1}^{lk} \mathbf{D}\_i e^{-\tilde{\gamma}\_i L} e^{-s\tau\_i} k = 1, 2, \dots, \text{Ng} \tag{34}$$

with *Ik* being the number of modal terms in the *k–th* group. To determine whether a set of exponential factors can be grouped or not, the maximum phase shifts associated to their time delays are compared. The set is grouped into a single delay group if the phase shift differences are below a pre-established value typically chosen at 10 (Morched et al., 1999). Each term **H***k* at (34) can now be considered free of delay factors and can be fitted as follows:

$$\tilde{\mathbf{H}}\_k = \sum\_{i=1}^{Nh(k)} \frac{\mathbf{R}\_{k,i}}{s - p\_{k,i}} k = 1, 2, \dots, Ng \tag{35}$$

where *Nh(k)* is the fitting order for the *k–th* term **H***k* , *pk,i* represents its *i–th* fitting pole and *Rk,i* is the corresponding matrix of residues. Note in (35) that common poles are being used to fit all elements at each matrix **H***k* . As (35) is introduced in (33), the following rational form is obtained for *H*:

$$\mathbf{H} = \sum\_{k=1}^{N\_{\mathcal{S}}} e^{-s\tau\_k} \sum\_{i=1}^{Nh(k)} \frac{\mathbf{R}\_{k,i}}{s - p\_{k,i}} \tag{36}$$

Initial estimates for the poles as well as time delays are obtained in the modal domain. The poles result from applying VF to the sum of the modal exponential factors conforming each delay group. The time delays proceed from Bode's formula as it has been said before. With all the poles *pk,i* and group delays *k* of (36) being estimated in the modal domain, the overall fitting of *H* is completed in the phase domain by obtaining the matrices of residues *Rk,i* and recalculating the poles (Gustavsen & Nordstrom, 2008). The fitting parameters so obtained can be taken as final or can be further refined by an iterative process. VF is applied throughout all these fitting tasks and detailed descriptions of these processes can be found in (Gustavsen & Nordstrom, 2008; Gustavsen & Semlyen, 1999).

#### **4.3. State-space analysis**

With the rational representation of *YC* and *H* state–space forms to evaluate *ish,0* and *iaux,0* arise naturally. Consider first the case of *ish,0*. Taking (22) and introducing (27) yields

$$\mathbf{I}\_{sh,0} = \mathbf{G}\_0 \mathbf{V}\_0 + \sum\_{i=1}^{Ny} \mathbf{W}\_i \tag{37}$$

where

276 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

*L* is a minimum phase–shift function (Bode, 1945) and

1

The time delays in (32) can be initially estimated by applying Bode's relation for minimum

provides the desired separation of *H* as a sum of terms, each one involving a single delay factor, the following consideration is brought in for computational efficiency (Morched et al., 1999). Modal delays often occur in groups with almost identical values. Suppose that a

1

*k*

*Ik L s*

( ) , 1 ,

*i k i*

*e*

where *Nh(k)* is the fitting order for the *k–th* term **H***k* , *pk,i* represents its *i–th* fitting pole and *Rk,i* is the corresponding matrix of residues. Note in (35) that common poles are being used to fit all elements at each matrix **H***k* . As (35) is introduced in (33), the following rational

 *s p* **R**

*Nh k k i*

with *Ik* being the number of modal terms in the *k–th* group. To determine whether a set of exponential factors can be grouped or not, the maximum phase shifts associated to their time delays are compared. The set is grouped into a single delay group if the phase shift differences are below a pre-established value typically chosen at 10 (Morched et al., 1999). Each term **H***k* at (34) can now be considered free of delay factors and can be fitted as follows:

*Ng s k*

*e*

1,2,..., *i i*

*e e k Ng*

1,2,...,

*k Ng*

( ) ,

*s p*

**R**

1 1 , *k Ng Nh k s k i k i k i*

number *Ng* of these groups can be formed and that (32) can be modified as follows:

1

*k i i*

*k*

phase complex functions (Bode, 1945) to the magnitudes of *exp(-*

factors at (30) can be further decomposed as follows:

associated to the velocity of the *i–th* mode. Hence:

where exp *i*

where *Ng* is less or equal to *N*, and

by comparing (33) and (32):

form is obtained for *H*:

1

where *Di* is the rank–1 matrix obtained from pre–multiplying the *i–th* column of *M* by the *i– th* row of *M–1*. Matrix *Di* in fact is an idempotent (Marcano & Marti, 1997). The exponential

; 1,2, , *i ii L Ls e ee i N*

*i N L i i*

*i i N L s i i*

*k*

*e e*

**H D** (30)

(31)

**H D** (32)

**H H** (33)

*k* is the representative delay for the *k–th* group. Clearly,

**H D** (34)

**H**(35)

**H** (36)

*i* is the time delay

*iL)* in (30). Although (32)

*e*

$$\mathbf{W}\_{i} = \frac{\mathbf{G}\_{i}}{\mathbf{s} - \mathbf{q}\_{i}} \mathbf{V}\_{0} \text{ } i = \mathbf{1} \text{ } \mathbf{2} \text{ } \dots \text{N} \text{y} \tag{38}$$

Application of the Inverse Laplace Transform to (37) and (38) gives the following continuous-time-state-space (CTSS) form for *ish,0*:

$$\mathbf{i}\_{sh,0} = \mathbf{G}\_0 \mathbf{v}\_0 + \sum\_{i=1}^{Ny} \mathbf{w}\_i \tag{39}$$

and

$$\frac{d\mathbf{w}\_i}{dt} = q\_i \mathbf{w}\_i + \mathbf{G}\_i \mathbf{v}\_0 \text{ \(i = 1, 2, \dots, Ny\)}\tag{40}$$

The CTSS form to evaluate *iaux,0*, is derived now. On replacing the fitted form of *H* given by (36) into (22):

$$\mathbf{I}\_{\text{aux},0} = \sum\_{k=1}^{\text{Ng}} \sum\_{i=1}^{\text{Nh}(k)} \mathbf{X}\_{k,i} \tag{41}$$

with

$$\mathbf{X}\_{k,i} = \frac{R\_{k,i}}{s - p\_{k,i}} I\_{r\emptyset,\iota} e^{-s\tau\_k} ; \quad \begin{array}{l} i = 1, 2, \ldots, N\text{y} \\ k = 1, 2, \ldots, N\text{g} \end{array} \tag{42}$$

Application of the Inverse Laplace Transform to (41) and (42) renders the following CTSS form:

$$\mathbf{i}\_{aux,0} = \sum\_{k=1}^{\text{Ng}} \sum\_{i=1}^{\text{Nh}(k)} \mathbf{x}\_{k,i} \tag{43}$$

$$\frac{d\mathbf{x}\_{k,i}}{dt} = p\_{k,i}\mathbf{x}\_{k,i} + \mathbf{R}\_{k,i}\mathbf{i}\_{rfl,L}(t-\tau\_k);\ \begin{array}{c} i = 1,2,...,Ny\\k = 1,2,...,Ny \end{array} \tag{44}$$

CTSS forms (39), (40), (43) and (44) provide the basis for a phase domain line model (Morched et al., 1999). Nevertheless, their solution by a digital processor requires the conversion to discrete–time state–space (DTSS). This is accomplished by applying a numerical differentiation rule to the CTSS forms. The one adopted here is the mid–point rule of differentiation, which is equivalent to the trapezoidal integration rule extensively used in EMTP (Dommel, 1969, 1992). Application of this rule to (44) with *t* as the solution time step results in:

$$\mathbf{x}\_{k,i} = a\_{k,i}\mathbf{x}\_{k,i}^{\prime} + \tilde{\mathbf{R}}\_{k,i}[\mathbf{i}\_{r\emptyset,L}(t-\tau\_k) + \mathbf{i}\_{r\emptyset,L}^{\prime}(t-\tau\_k)]; \quad \begin{array}{l} i = 1,2,...,Ng \\ k = 1,2,...,Ng \end{array} \tag{45}$$

where *ak,i*=(*2+tpk,i*)/ (*2-tpk,i*) and **R***k i*, =(*tRk,i*)/ (*2-tpk,i*). *xk,i* are discrete-state variables and primed variables denote their value at one previous time step *x*'*k,i= xk,i*(*t-t*). The discrete– time version of (43) maintains its original form:

$$\mathbf{i}\_{\text{aux},0} = \sum\_{k=1}^{N\_{\text{S}}} \sum\_{i=1}^{Nh(k)} \mathbf{x}\_{k,i}$$

Transmission line simulation of EMTs requires the use of time steps *t* smaller than any of the travel times *k* in the line. Hence, (45) provides the update of state vectors *xk,i* using only past values of variables already available, either from initial conditions or from previous simulation time steps.

The differentiation mid–point rule is now applied to (40):

$$\mathbf{w}\_{i} = a\_{i}\mathbf{w}\_{\ i}^{\prime} + \tilde{\mathbf{G}}\_{i}(\upsilon\_{0} + \upsilon\_{0}^{\prime}); i = 1, 2, \dots, N y \tag{46}$$

where *ai*=(*2+tqi*)/ (*2-tqi*) and **G***i* =(*tGi*)/ (*2-tqi*)

Expression (46) is not a proper DTSS form, as *wi* depends on the present–time value of *v0* which still is to be determined (Gustavsen & Mahseredjian, 2007). This problem is fixed here with the following redefinition of the state variable vector:

$$\mathbf{y}\_{i} = (\tilde{\mathbf{G}}\_{i}^{-1} \mathbf{w}\_{i} - \upsilon\_{0}) / (a\_{i} + 1); i = 1, 2, \dots, Ny \tag{47}$$

Introducing (47) in (46) and (39) the following expressions are obtained:

$$\mathbf{y}\_{i} = a\_{i}\mathbf{y}\_{i}^{\prime} + \boldsymbol{\upsilon}\_{\;\;0}^{\prime}; i = 1, 2, \ldots, N \\ y \tag{48}$$

$$\mathbf{\dot{i}}\_{sh,0} = \mathbf{G} \mathbf{v}\_0 + \mathbf{\dot{i}}\_{y-aux,0} \tag{49}$$

where

278 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

,, , ,

, ,, , , ,

primed variables denote their value at one previous time step *x*'*k,i= xk,i*(*t-*

Transmission line simulation of EMTs requires the use of time steps

The differentiation mid–point rule is now applied to (40):

*tqi*) and **G***i*

with the following redefinition of the state variable vector:

 =(*tGi*)/ (*2-*

1

Introducing (47) in (46) and (39) the following expressions are obtained:

*tpk,i*) and **R***k i*, =(

*k i k i k i rfl L k d i Ny p t dt k Ng*

CTSS forms (39), (40), (43) and (44) provide the basis for a phase domain line model (Morched et al., 1999). Nevertheless, their solution by a digital processor requires the conversion to discrete–time state–space (DTSS). This is accomplished by applying a numerical differentiation rule to the CTSS forms. The one adopted here is the mid–point rule of differentiation, which is equivalent to the trapezoidal integration rule extensively used in EMTP (Dommel, 1969, 1992).

> 1,2,..., ' [ ( ) ' ( )]; 1,2,..., *k i k i k i k i rfl L k rfl L k i Ny a tt k Ng*

> > *tRk,i*)/ (*2-*

past values of variables already available, either from initial conditions or from previous

0 0 ' ( ' ); 1,2,..., *i ii i* **w wG** *a v v i Ny* (46)

*tqi*)

Expression (46) is not a proper DTSS form, as *wi* depends on the present–time value of *v0* which still is to be determined (Gustavsen & Mahseredjian, 2007). This problem is fixed here

,

*k i*

**x**

Application of this rule to (44) with

*tpk,i*)/ (*2-*

time version of (43) maintains its original form:

where *ak,i*=(*2+*

the travel times

where *ai*=(*2+*

simulation time steps.

*tqi*)/ (*2-*

( ) ,0 , 1 1 *Ng Nh k aux k i k i*

1,2,..., ( ); 1,2,...,

*t* as the solution time step results in:

( ) ,0 , 1 1 *Ng Nh k aux k i k i* **i x**

*k* in the line. Hence, (45) provides the update of state vectors *xk,i* using only

0 ( ) / ( 1); 1,2,..., *i ii i v a i Ny* **y Gw** (47)

0 ' ' ; 1,2,..., *i ii* **y y** *a v i Ny* (48)

*sh*,0 0 ,0 *y aux* **i Gv i** (49)

 **x x Ri i** (45)

**i x** (43)

**x Ri** (44)

*tpk,i*). *xk,i* are discrete-state variables and

*t*). The discrete–

*t* smaller than any of

$$\mathbf{i}\_{y-aux,0} = \sum\_{i=1}^{Ny} \hat{\mathbf{G}}\_i \mathbf{y}\_i \; ; \; \hat{\mathbf{G}}\_i = (a\_i + 1)\tilde{\mathbf{G}}\_i \; ; \; \mathbf{G} = \mathbf{G}\_0 + \sum\_{i=1}^{Ny} \tilde{\mathbf{G}}\_i \mathbf{e}\_i$$

Expression (48) is a proper DTSS form for the sequential evaluation of *ish,0* at the phase domain line model.

Finally, the introduction of (43) and (49) in (23) results in:

$$
\dot{\mathbf{i}}\_0 = \mathbf{G} \mathbf{v}\_0 + \dot{\mathbf{i}}\_{hist, 0} \tag{50}
$$

with

$$\mathbf{i}\_{hist,0} = \mathbf{i}\_{y-aux,0} - \mathbf{i}\_{aux,0} = \sum\_{i=1}^{Ny} \hat{\mathbf{G}}\_i \mathbf{y}\_i - \sum\_{k=1}^{Ng} \sum\_{i=1}^{Nh(k)} \mathbf{x}\_{k,i}$$

Expression (50), along with (48) and (49), provides a discrete time–domain model for end *0* of the line segment at figure 1. The expressions for the model at end *L* are simply obtained by exchanging sub–indexes *0* and *L* at (48), (49) and (50). Obviously, state variables "*yi*" and "*xk,i*" of end *L* model are different from those of end *0*. Figure 3 provides a discrete–time circuit–model for the line segment of length *x=L*. This model is based on expression (50) and its companion for line end *L*. Note that the model consists of parallel arrangements of shunt conductances and auxiliary sources of currents comprising historic terms of ends (or nodes) *0* and *L*. Figure 3 model is thus in an appropriate form for computer code implementation. In this chapter, the Matlab environment has been chosen for this end.

**Figure 3.** Discrete time domain circuit representation of a multi-conductor line.

#### **5. Line model implementation in Matlab**

The discrete–time line model depicted in figure 3 and defined by (50) has been programmed by these authors in Matlab as an M–code function (see Appendix). This function consists of two sub–blocks, one for each multi-conductor line end. This model is to be used with a nodal network solver, a complete explanation on the nodal solver can be found in (Dommel, 1969 & 1992). Expression (50) constitutes essentially the interface between the line model and the nodal solver. Each one of the two sub–blocks in the line model performs iteratively the six tasks that are described next for line–end *0* sub–block. Figure 4 provides the block diagram of the complete line/cable model, along with its interfacing with the nodal solver.

$$\mathbf{i}\_{sh,0} = \mathbf{G}\mathbf{v}\_0 + \mathbf{i}\_{y-au\alpha,0}$$

**Step 3.** Auxiliary current source value, due to the reflected traveling waves at the remote line end, is updated by (43):

$$\mathbf{i}\_{\text{aux},0} = \sum\_{k=1}^{N\_{\text{S}}} \sum\_{i=1}^{Nh(k)} \mathbf{x}\_{k,i}$$

**Step 4.** Vector of reflected currents at the local line–end (node) "*irfl,0*" is calculated for the present time by means of (26) being modified to suit line–end *0*:

$$\mathbf{i}\_{r\emptyset,0} = \mathbf{2i}\_{sh,0} - \mathbf{i}\_{au\alpha,0}$$

This vector is delivered to end *L* sub–block through a delay buffer. Although branch current vector *i0* usually is not explicitly required, it is conveniently evaluated here by (50):

$$\mathbf{i}\_0 = \mathbf{G}\mathbf{v}\_0 + \mathbf{i}\_{hist, 0}$$

**Step 5.** Internal states inside the line model are updated by (48) and (45):

$$\mathbf{y}\_i = a\_i \mathbf{y}'\_i + \boldsymbol{\upsilon}'\_0$$

$$\mathbf{x}\_{k,i} = a\_{k,i} \mathbf{x}'\_{k,i} + \tilde{\mathbf{R}}\_{k,i} [\mathbf{i}\_{r\emptyset,L}(t - \boldsymbol{\tau}\_k) + \mathbf{i}'\_{r\emptyset,L}(t - \boldsymbol{\tau}\_k)]$$

**Step 6.** The vector of history currents for end (node) *0* is updated by means of (50) and the update is delivered to the nodal-network solver.

Steps *1* to *6* are iterated *Nt* times until *Ntt* spans the total simulation time of interest.

#### **5.1. Handling of line-travel delays**

It follows from expressions (43) and (45) that the calculation of *iaux,0* requires the reflected currents vector *irfl,L* being evaluated with various time delays , …, *Ng*. Recall that the delays are due to the travel time needed by a wave to travel from one line end to the other. Past values of *irfl,L* can be obtained either from line initial conditions or from previous simulation steps; nevertheless, these values are given by samples regularly distributed *t* seconds apart. Since the involved travel times (or line delays) usually are not integer multiples of *t*, the required values of *irfl,L* must be obtained by means of interpolations. The standard procedure for this is to interpolate linearly (Dommel, 1992) and this is adopted here.

Evaluation of the delayed values require a memory buffer spanning at least the largest travel time

$$
\pi\_{\text{max}} = \max \{ \tau\_{1'} \tau\_{2'} ... \tau\_{\text{Ng}} \} \, \, \, \tag{51}
$$

and buffer length *Nb* is calculated as follows:

280 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 1

to determine line end (nodal) voltages *v0* and *vL*.

repeated here for convenience:

line end, is updated by (43):

**Step 1.** State–variable and history–current values are assumed known, either from initial conditions or from previous simulation steps. These values are used by the nodal solver

**Step 2.** Shunt current due to the characteristic admittance of the line is calculated by (49)

*sh*,0 0 ,0 *y aux* **i Gv i**

**Step 3.** Auxiliary current source value, due to the reflected traveling waves at the remote

**Step 4.** Vector of reflected currents at the local line–end (node) "*irfl,0*" is calculated for the

,0 ,0 ,0 2 *rfl sh aux* **i ii**

This vector is delivered to end *L* sub–block through a delay buffer. Although branch current

0 0 ,0 *hist* **i Gv i**

0 ' ' *i ii* **y y** *a v*

, ,, , , , ' [ ( ) ' ( )] *k i k i k i k i rfl L k rfl L k* **x x Ri i** *a tt*

**Step 6.** The vector of history currents for end (node) *0* is updated by means of (50) and the

It follows from expressions (43) and (45) that the calculation of *iaux,0* requires the reflected

delays are due to the travel time needed by a wave to travel from one line end to the other. Past values of *irfl,L* can be obtained either from line initial conditions or from previous simulation steps; nevertheless, these values are given by samples regularly distributed

seconds apart. Since the involved travel times (or line delays) usually are not integer

standard procedure for this is to interpolate linearly (Dommel, 1992) and this is adopted

*t*, the required values of *irfl,L* must be obtained by means of interpolations. The

*t* spans the total simulation time of interest.

, …,

*Ng*. Recall that the

*t*

vector *i0* usually is not explicitly required, it is conveniently evaluated here by (50):

present time by means of (26) being modified to suit line–end *0*:

**Step 5.** Internal states inside the line model are updated by (48) and (45):

update is delivered to the nodal-network solver.

currents vector *irfl,L* being evaluated with various time delays

Steps *1* to *6* are iterated *Nt* times until *Nt*

**5.1. Handling of line-travel delays**

multiples of

here.

( ) ,0 , 1 1 *Ng Nh k aux k i k i* **i x**

$$N\_b = \left\lfloor \frac{\tau\_{\text{max}}}{\Delta t} \right\rfloor + 1\tag{52}$$

**Figure 4.** Line/Cable model's complete flow diagram.

If a propagation delay is an integer multiple of *t*, the required value of *irfl* can be readily retrieved from the memory buffer. This is illustrated by figure 5 where the simulation time step is *t=0.03* ms and the travel time is =0.10 ms. It can be seen that at simulation time *t= 0.24* ms the required history value at 0.09ms is available from the table.

On a multiphase system, nevertheless, it is highly improbable that all the propagation times can be made integer multiples of a single value of *t* suitable for transient simulations. Thus, the required values must be obtained through interpolation. Figure 6 illustrates this case, where a simulation time step *t=0.04* ms is assumed instead of the *t=0.03* ms one at figure 5. Notice that now the required history value, for a time delay of 0.09ms, is not readily available.

Suppose now that the required value *irfl(t–)* is between the *k–th* and the *(k+1)–th* stored samples of *irfl*. Let be the fractional part of */t*, that can be obtained as follows:

$$
\zeta = \frac{\pi}{\Delta t} - \left\lfloor \frac{\pi}{\Delta t} \right\rfloor,\tag{53}
$$

with, 0 < <1. The estimated value of *irfl(t–)* by linear interpolation is thus:

$$\mathbf{\dot{i}}\_{r\|,L}(t-\tau\_k) \equiv \mathbf{\dot{i}}\_{r\|,L}(t-r\Delta t) + \boldsymbol{\zeta}[\mathbf{\dot{i}}\_{r\|,L}(t-k\Delta t) - \mathbf{\dot{i}}\_{r\|,L}(t-k\Delta t-\Delta t)].\tag{54}$$

Figure 7 illustrates the memory buffer management, either for *irfl,0* or for *irfl,L*. At the first simulation time step corresponding to time *t=0t*, calculated *irfl* is stored at memory *1*, and so on until step *Nb* which is the buffer size limit. Beyond this limit, memory cells *1*, *2*, *3* and on, are overwritten as figure. 7 shows, since their previously stored values are not needed any longer.

**Figure 5.** Interpolation scheme: *t* integer multiple.

**Figure 6.** Interpolation scheme: *t* non integer multiple.

**Figure 7.** History buffer management.

be the fractional part of

where a simulation time step

Suppose now that the required value *irfl(t–*

<1. The estimated value of *irfl(t–*

simulation time step corresponding to time *t=0*

available.

with, 0 <

any longer.

**Figure 5.** Interpolation scheme:

**Figure 6.** Interpolation scheme:

*t* integer multiple.

*t* non integer multiple.

samples of *irfl*. Let

the required values must be obtained through interpolation. Figure 6 illustrates this case,

5. Notice that now the required history value, for a time delay of 0.09ms, is not readily

*/*

,, , , ( ) ( ) [ ( ) ( )]. *rfl L k rfl L rfl L rfl L* **ii i i** *t t rt t kt t kt t*

Figure 7 illustrates the memory buffer management, either for *irfl,0* or for *irfl,L*. At the first

so on until step *Nb* which is the buffer size limit. Beyond this limit, memory cells *1*, *2*, *3* and on, are overwritten as figure. 7 shows, since their previously stored values are not needed

> *t(ms) irfl,0 irfl,L* 0.0 \* \* 0.03 \* \* 0.09 \* \* 0.12 \* \* 0.15 \* \* 0.18 \* \* 0.21 \* \* 0.24 \* \*

*t(ms) irfl,0 irfl,L* 0.0 \* \* 0.04 \* \* 0.08 \* \* 0.12 \* \* 0.16 \* \* 0.20 \* \* 0.24 \* \*

*t=0.04* ms is assumed instead of the

, *t t*

 *)* is between the *k–th* and the *(k+1)–th* stored

*t*, calculated *irfl* is stored at memory *1*, and

*t*, that can be obtained as follows:

*)* by linear interpolation is thus:

(53)

*t=0.03* ms one at figure

(54)

Linear interpolation is an order 1 numerical procedure and the trapezoidal rule used for the rest of the line model is of order 2. The question arises as to whether or not the order 2 quadratic interpolation should be adopted instead. This has been investigated at (Gutierrez– Robles et al., 2011) and it has been found that the increase in accuracy is marginal.


*Editorial* **Advances in Mineral Processing and Hydrometallurgy**

**Corby G. Anderson \* and Hao Cui \***

Kroll Institute for Extractive Metallurgy, Mining Engineering Department & George S. Ansell Department of Metallurgical and Materials Engineering, Colorado School of Mines, Golden, CO 80401, USA **\*** Correspondence<EMAIL_ADDRESS>(C.G.A<EMAIL_ADDRESS>(H.C.)