**3. Biology**

*Goats (Capra) - From Ancient to Modern*

data of animal numbers.

enough attention to collect exact data.

• Cashmere is considered more likely as a by-product, so governments do not pay

• The production of cashmere cannot be estimated from the number of goats. There is no direct relation between the number of animal or animal weight and the production of cashmere. For the reason of existing very different cashmere-producing breeds/ecotypes with diverse genetic potential of cashmere production, it is difficult to estimate total cashmere production from the

*Goats and sheep population and cashmere prices in Mongolia, 1992–2000, Sources: Bakey et al. (2001); FAO (2002).*

*Raw cashmere purchase price from farmers (LPP, LIFE Network, IUCN–WISP and FAO. 2010. Adding value to livestock diversity – Marketing to promote local breeds and improve livelihoods. FAO Animal Production*

**86**

**Figure 3.**

**Figure 2.**

*and Health Paper. No. 168. Rome).*


**1. Introduction**

Since the beginning of human civilization, construction and demolition wastes (CDW) had been becoming a global problem that affects the sustainable development of the resources and the environment [1]. To date, more than a hundred billion tons of CDW were generated in the world, while approximately 30% to 50% of them were waste concrete [2]. However, this considerable amount of waste concrete is mainly disposed of in landfills, resulting in a severe environmental problem [3–6]. In order to protect the ecological environment and achieve sustainable development, reducing, reusing, and recycling construction waste is a desirable approach to preserve the ecological environment [7–11]. Therefore, many countries have passed legislation to encourage the recycling of waste concrete for the resource utilization of recycled concrete.

Recycled concrete aggregates (RCA) can be obtained from waste concrete by crushing, screening, cleaning, and separating [12,13]. Previous studies have been conducted to characterize the potential advantages and drawbacks of RCA. Being different from natural aggregates (NA), RCA has a layer of old mortar attached to the surface, which has more loose pores at the interface [14]. The attached old mortar brings with it worse properties: a lower apparent density and higher water absorption [15,16]. Consequently, these properties cause the mechanical properties of RCA to be inferior to that of NA, which limits the application of recycled concrete in civil engineering. Although the mechanical properties of recycled aggregates are poor, different scholars, institutions, and countries have been applying recycled aggregates into structures through rational design and experiment in recent years, and in some cases, good results have been achieved [17–23].

To further promote the recycling of construction waste, growing studies on the production of recycled concrete blocks prepared by RCA have been conducted in years. Researchers proposed that the mechanical properties of recycled concrete blocks mainly depend on the substitution rate and qualities of RCA. Soutsos et al. [24] pointed out that the increase of recycled fine aggregates content has a more significant impact on the reduction of recycled concrete block strength. Therefore, the maximum replacement rate of recycled fine aggregates is recommended to be 20%. Bai et al. [25] also confirmed the findings through a similar experiment. A study conducted by Poon et al. [26] aimed to investigate the mechanical properties of recycled concrete blocks. They found that the substitution rate of RCA below 50% had little e ffect on the compressive strength of recycled concrete blocks. However, the compressive strength of recycled concrete hollow block (RCHB) decreases with the increase of the content of low-grade RCA (The content of soil or broken brick in aggregates being >10%) [27]. Guo et al. [28] investigated the influence of di fferent substitution rates of recycled concrete aggregates on the mechanical properties of recycled concrete blocks. The test results illustrated that with the substitution rate up to 75%, the strength of recycled concrete blocks slightly decreases but still complies with the standards. An experiment by Sabai et al. [29] showed that the compressive strength of recycled concrete blocks with 100% RCA can achieve the target of 7 MPa, and even the minimum strength requirement of construction.

Despite the fact that the compressive strength of the concrete block is the critical performance index about whether the block can be used in a masonry structure, the mechanical performance of a masonry prism is more reflective of the actual stress state of the masonry structure. Corinaldesi et al. [30] studied the compressive, shear, and bond strength of recycled mortar prisms, and found that the shear and compressive strength of recycled mortar prisms were close to or even better than that of ordinary mortar prisms. Guo et al. [28] conducted a study of recycled concrete block prisms, and the conclusions are consistent with that of Corinaldesi et al. [30]. In general, the shear strength of the recycled concrete prisms is close to that of the ordinary concrete prisms. However, the seismic performance of the masonry structure, especially the hysteretic characteristics and energy dissipation capacities of the masonry is unclear by the tests of the masonry prisms.

Through the analysis of the above examples, it can be proved that the mechanical properties of recycled blocks can satisfy the requirement of practical application, and the shear properties of masonry assemblages fabricated by recycled concrete blocks are close to those of ordinary masonry assemblages. However, there is still a lack of knowledge about the recycled concrete blocks, and whether they can be used to produce masonry structures with seismic requirements, especially regarding the research about the seismic performance of masonry structures, which is crucial for addressing the utilization problem of RCA. Therefore, if the waste concrete can be recycled to produce RCHB which can be used in the structures with seismic requirements, this type of RCHB will be popularized and further applied on a larger scale.

In view of this, RCA generated from waste concrete was used to produce a new type of RCHB which can be used for masonry structures. In this experiment, three RCHB masonry walls without structural columns constraint, one RCHB wall constrained by the recycled aggregate concrete structural columns, and other walls constrained by ordinary concrete structural columns were manufactured for seismic performance testing. The seismic behavior of the specimens, such as the failure pattern, the hysteresis curves, the skeleton curves, the ductility coe fficient, and the energy dissipation of the specimens was analyzed under cyclic loading. The influences of aspect ratio, vertical axial stress, and di fferent materials used for structural columns on the seismic performance of RCHB masonry walls were also studied. Finally, the seismic capacity of RCHB masonry structure under seismic loading was attained.


**Albana Kona 1,\*, Paolo Bertoldi 1 and ¸Siir Kılkı¸s 2**

Received: 30 December 2018; Accepted: 3 March 2019; Published: 13 March 2019

**Abstract:** Local authorities and cities are at the forefront of driving the energy transition, which plays a crucial role in mitigating the effects of climate change. The greenhouse gas emissions in cities, due to energy consumption, are placed into two categories: direct emissions generated from the combustion of fossil fuels mainly in buildings and transport sectors, and indirect emissions from grid-supplied energy, such as electricity and district heating and/or cooling. While there is extensive literature focused on direct greenhouse gas emissions accounting in cities' inventories, research has focused to a lesser extent on allocation methods of indirect emissions from grid-supplied energy. The present paper provides an updated definition for the concept of local energy generation within the Covenant of Mayors initiative and proposes a new methodology for indirect emission accounting in cities' greenhouse gas emission inventories. In addition, a broader policy framework in which local action is taken is discussed based on the European Union energy and climate policies, and over 80 exemplary Covenant of Mayors good practices are identified across the technology areas of local energy generation and four modes of urban climate governance. The contributions of the paper demonstrate that local authorities have the capacity to support and mobilize action for local energy generation investments through the multiple modes of urban climate governance to update and strengthen climate action

**Keywords:** Covenant of Mayors; indirect emission allocation; local energy generation


**Author details**

Agnes Chinecherem Nkele1 \* and Fabian I. Ezema1,2,3,4

1 Department of Physics and Astronomy, University of Nigeria, Nsukka, Enugu, Nigeria

2 Nanosciences African Network (NANOAFNET), iThemba LABS-National Research Foundation, 1 Old Faure Road, Somerset West 7129, P.O. Box 722, Somerset West, Western Cape Province, South Africa

3 UNESCO-UNISA Africa Chair in Nanosciences/Nanotechnology, College of Graduate Studies, University of South Africa (UNISA), Muckleneuk Ridge, P.O. Box 392, Pretoria, South Africa

4 Department of Physics, Faculty of Natural and Applied Sciences, Coal City University, Enugu, Nigeria

\*Address all correspondence to<EMAIL_ADDRESS>
© 2020 The Author(s). Licensee IntechOpen. 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.


**4. Targeting HIF-responsive growth factors**

## **4.1. Vascular endothelial growth factor**

Clear cell renal carcinomas are notoriously angiogenic. Indeed, prior to the availability of computed tomography, renal angiograms were often used to diagnose these tumors. Renal carcinomas overproduce a variety of angiogenic moieties including vascular endothelial growth factor (VEGF), the product of a HIF-responsive gene. In addition to promoting angiogenesis, VEGF might suppress antitumor immune responses as well. It has also been suggested that VEGF has direct stimulatory effects on renal carcinoma cells, although these findings await further corroboration, [8, 9].

Several drugs that inhibit VEGF, or its kinase insert domain-containing receptor (KDR), have activity against clear cell renal carcinomas. In a randomized phase II study, patients with metastatic renal carcinoma who were treated with 10 mg/kg (but not 3 mg/kg) bevacizumab, a neutralizing antibody against VEGF, exhibited a significant delay in time-to-disease pro‐ gression, [10]. Other unrelated KDR inhibitors such as SU11248 (sunitinib maleate), BAY43-9006 (sorafenib), and AG-013736 also appear to have significant activity against this tumor subtype, [11, 12].


**2.1 Physical therapy**

Physical therapy relieves pain and muscle stiffness, improves mobility, and prevents future problems such as muscle tightness. It is always recommended as the first therapy intervention. It improves motor skills and can stop movement problems from worsening in the future. It includes varieties of treatment modalities such as strength, flexibility exercises, heat treatment, massages, and special equipment for CP children to give them more independence. How much the PT will help to improve the child's condition depends on the severity and CP type the child has. Mild cases only require some physical therapy. The best chance to improve the child's condition is with early intervention. Mobility can be improved by preventing contractures and joint dislocations by keeping the body strong and flexible. PT can improve the level of self-reliance of most CP children. The goal of PT is to make everyday movements easier. It can improve coordination, balance, strength, flexibility, endurance, pain management, posture, gait, and overall health.

**The benefits of PT depend on the CP type:**

The treatment plan will be based on the CP type. The child might have either hemiplegia, diplegia, or quadriplegia. Therapists design specific exercises and routines for each type, which helps the child regain movement in the affected area over time.

Other conditions such as scoliosis, thoracic kyphosis, lumbar lordosis, pelvic inclination, pelvic rotation, pelvic obliquity, knee deformity, shortened Achilles tendon, and hand and wrist deformities can be treated by PT. After the initial evaluation by the therapist, a treatment plan will be created, which might include a combination of exercises, muscle relaxing techniques, and special equipment to improve movement.

CP child either will have muscle hypotonia or hypertonia. Flexibility exercises and massages are usually used for the hypertonic child, which will develop mobility and prevent painful muscle tightening, which could require surgical intervention. On the other hand, the athetoid child will require strength training exercises to increase their muscle tone. Walking, posture, transitional movements, and sensory impairments such as touch and balance can be improved by certain exercises. Exercises that improve sitting, kneeling, and standing can improve posture. Infants use transitional movements to be able to walk, roll over, and sit up are examples of those movements.

*Management of Spasticity and Cerebral Palsy Update DOI: http://dx.doi.org/10.5772/intechopen.106351*

Different mobility aids are to be used to make therapy more effective. Different orthotics equipment can be used to help with walking, posture, and joint mobility; braces, casts, splints, and shoe inserts are some examples. There are other tools that have been used in PT include exercise balls, resistance bands, free weights, swimming pools, hot and cold packs, and electric muscle stimulation (it improves gait and upper limb function).

As the child gets older, the physical therapist will change the program to meet the child's needs. Toddlers and the early school years are the most important stages of a child's life.

Outcome of the therapy can be evaluated by observation and a set of standardized tests to measure motor function. Therapists will be looking for specific postures and movements that could be corrected to develop a strong treatment plan.

#### **2.2 Occupational therapy**

Health care providers usually recommend OT along with PT for CP children to provide easier independent lives because it has many benefits. On one hand, PT improves muscle tone and balance, and on the other hand, OT is helping children to accomplish daily living activities. According to the American Occupational Therapy Association (AOTA), OT helps people of all ages to learn how to improve their abilities to do their everyday activities [3]. Everyday school tasks, eating and drinking, dressing, hygiene, and other tasks necessary for daily use can be assisted by OT. Children with visual impairment, regardless of severity, can do their daily tasks by doing OT. Also, sensoryperceptual-motor (SPM) training, which helps to stimulate the senses of performing daily tasks more easily, can be assisted by OT. An occupational therapist usually uses visual input tools and cues, as well as verbal directions to help children in carrying out necessary everyday tasks. A cognitively impaired child can benefit from OT. Establishing daily schedules or routines to help those children some daily tasks such as remembering to dress, brush their teeth, and other tasks are needed to carry out their day.

The real challenge for CP children is their motor impairment, which makes their daily tasks difficult. For example, most children have difficulties with daily hygiene, moving around, posture, and completing school tasks due to their movement limitations. The occupational therapist will assess every child's needs individually to design a program to meet their personal goals. The main line of OT program includes stretching and guide motor output techniques and constraint-induced movement therapy (CIMT), which helps hemiplegic CP children [4], and neuromuscular facilitation techniques.

OT can help in other areas such as grasping objects; using computer and other electronics; opening doors; writing, holding, and reading books; playing; interacting with peers, caregivers, and parents; school and/or daycare activities; and sitting and standing up from tables and school desks. CP Children and their parents and caregivers can benefit from OT by watching children thrive both physically and emotionally. This progress reduces parents' stress and provides a sense of peace as the child becomes more independent and productive.

OT can take place in different settings such as outpatient offices, inpatient medical centers, at home, nursing centers, and OT providers. To obtain quality OT, the child must be under the care and supervision of a trained and licensed occupational therapist. By evaluating the child and assessing the child's strengths and weaknesses, they will be able to formulate an individualized plan to help the child best reach his or her goals. Per AOTA, licensed occupational therapists, besides they are licensed in their state to practice occupational therapy, have another degree in Social Services, Health, Biology, Chemistry, Psychology, Sociology, Human Development, and Anatomy & Physiology. To obtain their license, they must be graduated from an accredited school, do fieldwork hours, and pass the National Board for Certification in Occupational Therapy examination. The therapist will conduct a comprehensive assessment including testing the child's: motor skills, developmental condition, cognitive development and functioning, psychological needs, emotional needs, and home and school environment to be able to create an individualized treatment plan.


*5.2.1. Styrene hydrogenation*

**Figure 10** shows a scheme of the reaction of selective hydrogenation of styrene to ethylbenzene. In this reaction the hydrogenation of the aromatic nucleus is undesired.

**Figure 10.** Scheme of the reaction of selective hydrogenation of styrene to ethylbenzene.

The reaction of hydrogenation of styrene to ethyl benzene is used as a reaction test for the process of purification of some refinery hydrocarbon streams. In this reaction the interest is to hydrogenate the vynillic external C═C bond, while keeping unaltered the bonds of the aromatic ring [34, 35].

In order to determine the mass transfer limitations in the inside of the catalyst particles, the effectiveness factor (η) was used:

$$\eta = \frac{\text{real reaction rate}}{\text{unstretrected reaction rate}} \tag{1}$$

The real reaction rate corresponds to the measured rate of the reaction at the conditions of temperature and concentration of the bulk reaction medium.The unrestricted reaction rate corresponds to the ideal reaction rate of a catalyst particle with no intraparticle or fluidparticle gradients of concentration or temperature that has the temperature and concentration of the bulk fluid phase.The unrestricted reaction rate was measured by grinding the catalyst pellets to a particle size smaller than 100 μm. For this particle size and the reaction conditions used it was estimated that all mass and heat transfer gradients were negligible.

The results of styrene conversion as a function of reaction time for the pelletized and powder catalysts are shown in **Figure 11**. The catalysts prepared with the composite supports UTAl and BTAl have higher activities (higher conversion values) than the catalysts prepared with commercial supports.

**Figure 11.** Styrene conversion as a function of reaction time.Catalysts with 1% Pd. Reaction conditions: toluene solvent, *C*0 styrene = 0.445 M, 2.0 MPa H2, 353 K, *W*cat = 0.3 g (pellet) or 0.01 g (powder).

In all the tests, for all catalysts, at all reaction conditions used, the obtained values of selectivity to ethylbenzene were higher than 99%.

Initial reaction rates were calculated by calculating the conversion rate as a function of time and extrapolating to zero conversion. With the values corresponding to the pellet and powder catalyst the effectiveness fact or could thus be estimated using the following equation:

New Strategies for Obtaining Inorganic-Organic Composite Catalysts for Selective Hydrogenation http://dx.doi.org/10.5772/65959 195

$$\eta = \frac{\text{Initial reaction rate of pellets catalyst}}{\text{Initial reaction rate of powder catalyst}} = \frac{r\_{p\text{ollder}}^0}{r\_{p\text{ouder}}^0} \tag{2}$$

**Table 3** contains the values of pellet 0 , powder 0 , and for the different catalysts. It can be seen that the composite supports have highervalues than PdRX and PdAl catalysts. The higher effec‐ tiveness values are directly related to the lower mass transfer limitations when using the composite catalysts.

hydrogenate the vynillic external C═C bond, while keeping unaltered the bonds of the

In order to determine the mass transfer limitations in the inside of the catalyst particles, the

The real reaction rate corresponds to the measured rate of the reaction at the conditions of temperature and concentration of the bulk reaction medium.The unrestricted reaction rate corresponds to the ideal reaction rate of a catalyst particle with no intraparticle or fluidparticle gradients of concentration or temperature that has the temperature and concentration of the bulk fluid phase.The unrestricted reaction rate was measured by grinding the catalyst pellets to a particle size smaller than 100 μm. For this particle size and the reaction conditions used it

The results of styrene conversion as a function of reaction time for the pelletized and powder catalysts are shown in **Figure 11**. The catalysts prepared with the composite supports UTAl and BTAl have higher activities (higher conversion values) than the catalysts prepared with

**Figure 11.** Styrene conversion as a function of reaction time.Catalysts with 1% Pd. Reaction conditions: toluene solvent,

In all the tests, for all catalysts, at all reaction conditions used, the obtained values of selectivity

Initial reaction rates were calculated by calculating the conversion rate as a function of time and extrapolating to zero conversion. With the values corresponding to the pellet and powder catalyst the effectiveness fact or could thus be estimated using the following equa-

styrene = 0.445 M, 2.0 MPa H2, 353 K, *W*cat = 0.3 g (pellet) or 0.01 g (powder).

to ethylbenzene were higher than 99%.

= (1)

real reaction rate unrestricted reaction rate

h

194 New Advances in Hydrogenation Processes - Fundamentals and Applications

was estimated that all mass and heat transfer gradients were negligible.

aromatic ring [34, 35].

commercial supports.

*C*0

tion:

effectiveness factor (η) was used:

**Table 3.** Styrene hydrogenation initial reaction rate for catalyst pellets (*r*0 pellet) and catalysts powders (r0 powder). Effectiveness factor (η) [31].

Catalysts with a metal loading lower than 0.5 Pd wt% were also tried with this reaction test. The catalysts are named 0.3PdUTAl and 0.3PdBTAl in the case of the composites (with 0.3% Pd). LD265 is a procatalyse Pd/Al2O3 catalyst (0.3% Pd) and Engelhard is a BASF Pd/Al2O3 catalyst (0.5% Pd).The results are plotted in **Figure 12**.

**Figure 12.** Styrene hydrogenation catalysts with less of 0.5 wt % of Pd. Reaction conditions: toluene solvent, *C*0 styrene = 0.445 M, 2.0 MPa H2, 353 K, *W*cat = 2 g pellet.

It can be seen that for all catalysts the pattern of conversion as a function of time are quite similar, though the catalytic activity is slightly better for the composite 0.3PdUTAl and 0.3PdBTAl catalysts.


**2.2 PON1 and dementias**

Erstwhile dementia is now termed as major neuro-cognitive disorder (MND). It describes an overall decline in memory and cognitive skills severe enough to reduce person's ability to perform everyday activities. Previously thought to affect elderly only now it is affecting the younger age group too. It is characterized by significant decline in any of the cognitive domain including, executive function, complex attention, language, learning, memory, and perceptual-motor or social cognition. Major neuro-cognitive disorder is diagnose by the decline in patient's previous cognitive ability without delirium, which should be persistent and progressive over the time. At present disease burden supposed to be 43 million worldwide which is expected to escalate to 131 million by 2050. Alzheimer disease (AD) which is responsible for nearly 70–80% cases of dementia worldwide, is one of the important cause of death over the age of 65 years [31].

One of the characteristic neuro-pathological features of AD is the presence of amyloid plaques which comprise aggregates of β-amyloid derived from the amyloid precursor protein [32]. Increasing evidence suggests that cholesterol plays a role in the pathophysiology of Alzheimer's disease, and elevated serum total-cholesterol level has been shown to be a risk factor for AD [33]. Abnormal phosphorylation of tau proteins is thought to be responsible for pathogenesis of AD. In addition to phosphorylation, neuronal degeneration is caused by a combination of beta amyloid production, oxygen deficiency. Lipid peroxidation too plays role in pathogenesis of AD. Amyloid is responsible for oxidative stress through free radicals. This oxidative stress is responsible for conversion of soluble amyloid proteins to insoluble fibrils and further in to polymerization of tau proteins [34, 35]. Some researchers report that deficient serum paraoxonase activity is a significant risk factor for AD and that paraoxonase activity is governed in part by at least 2 distinct variants, one located in the PON1 region and another in PON2 [36]. Some reports suggest that low PON1 activity is associated with cognitive decline, especially in AD [37]. There are many opposite results regarding polymorphism In MND. Some says there is definite correlation between the Q192R and L55M polymorphism and risk of AD and MND while others differ in their views as they found no correlation [38, 39].

### **2.3 PON1 and Parkinson's disease**

Parkinson's disease (PD) is an idiopathic disease of the nervous system characterized by both motor and non-motor system manifestations. It is a chronic progressive neurodegenerative disorder that occurs mostly in older persons but that can appear in much younger patients. It is the second most common neurodegenerative disease. Sometimes called "paralysis agitans", PD is uncommon in young people, especially those under 40 years of age [40, 41]. The pathological definition of PD is loss or degeneration of the dopaminergic neurons in the substantia nigra and development of Lewy Bodies in dopaminergic neurons. Pathologic changes may precede obvious symptoms by two decades or more. This preferential loss of dopamine producing neurons and simultaneous lack of cholinesterase inhibition results in marked imbalance acetylcholine and dopamine along with impairment of motor control. Lewy Bodies, or abnormal intracellular aggregates, contain various proteins including α-synuclein and ubiquitin that impair optimal neuron functioning [42, 43].

While the exact pathogenesis of the disease has not been completely elucidated, several theories implicating the association of genetic and environmental toxic elements such as exposure to pesticides or an oxidative cell environment conspire to trigger the neuron degeneration. Recent reports suggest that the chronic low-grade inflammation due to various sources like pesticides, drugs, aging process etc. is responsible for the cellular senescence in nervous tissue [44, 45]. From a pathologic

*Paraoxonase in Nervous System DOI: http://dx.doi.org/10.5772/intechopen.110843*

perspective, the brain's substantia nigra pars compacta and pontine locus ceruleus are affected by typical abnormalities of PD including depigmentation, neuronal loss and gliosis. Nearly 60–70% neurons are lost by the time symptoms appears [46, 47].

Pesticides have been implicated in the development of PD by inhibition of ubiquitine proteosome system. Probably upon exposure to pesticides, mitochondrial dysfunction and α-synuclein, which is a neuronal protein, undergo conformational change which leads to symptoms of PD. One hypothesis states that as a result of mitochondrial dysfunction, there is release of cytochrome-c in the cytosol which binds to apaf1 and starts apoptosis. This results into neuronal degeneration and ultimately to PD [48].

As PON1 has antioxidant capacity and is capable of hydrolyzing toxic substances, literature links PON1 polymorphism with PD. PON1 polymorphism as Met54 may be an independent risk factor for PD. This mutation could possibly be responsible for decreased PON1 activity which results in lessened metabolism of environmental neurotoxins and could play a role in neuro-degeneration. Researchers have found that there is an association between L55M polymorphism of PON1 and PD, whereas Q192R polymorphism was unlikely to be a major risk factor for susceptibility to PD [49]. Recently a researcher found MM PON1–55 genotype exhibit greater than 2-fold increase in PD risk when exposed to organophosphates, compared with subjects who had the wild type or heterozygous genotype and no exposure [50].

#### **2.4 PON1 and multiple sclerosis**

Multiple sclerosis (MS) is one of the most common neurological disorders, occurring mainly in young adults in age 20–40 years and more commonly in women. MS is a chronic inflammatory disease characterized by demyelinating lesions in the brain, spinal cord, and optic nerve [51, 52]. Some studies have suggested a role of the oxidative stress and lipid peroxidation in the pathogenesis of MS. Due to the pathogenetic role of reactive oxygen species and the oxidation of lipoproteins in MS pathology, antioxidants prevent free-radical mediated tissue destruction [53]. PON1 is shown to play an important antioxidant role in the blood. Variability of PON1 activity depends on polymorphism in the coding region. There are only a few studies that describe the relationship between PON1 polymorphism and the risk for MS [54]. The relationship between MS and oxidative stress due to oxidized lipoproteins in cholesterol transport has been discussed. It has been proposed that increased number of oxidized HDL particles due to decreased PON activity may be unable to protect LDL against oxidation increasing the risk of atherogenesis which may lead to MS [55].

Oxidative stress is a critical factor in pathogenesis of MS as it promotes leukocyte migration, participates in oligodendrocyte damage and axonal injury. Reactive oxygen species and reactive nitrogen species are produced in the CNS of MS patients chiefly by activated macrophages and microglia could account for demyelination and axonal disruption, the hallmarks of the disease [56, 57].

#### **2.5 PON1 and amyotrophic lateral sclerosis**

Amyotrophic lateral sclerosis (ALS) is a dangerous neurodegenerative disease, characterized by progressive motor neurons loss, paralysis, and inflammation with an average survival of 3–5years after diagnosis. 5–10% of cases are familial, and70% of the familial cases can be explained by identified gene mutations, e.g. the *C9orf72* repeat expansion [58]. ALS is a degenerative disease of adult-onset

and fatal outcome characterized by the simultaneous loss of motor neurons in the cerebral cortex, brainstem and spinal cord. In the familial cases, a large number of mutations have been found, most prominent of which is the mutation in the gene that encodes for the cytosolic enzyme superoxide dismutase [59, 60]. It was found in recent studies that PON1 gene is down-regulated in the central nervous system and peripheral cell like lymphocytes and fibroblasts. Some genetic mutations also are identified in the patients of ALS [61, 62]. It is well established fact that the oxidative stress is responsible for degeneration of the motor neurons in ALS. The CNS as a whole is particularly susceptible to oxidative stress because the neuronal membrane contains a high abundance of polyunsaturated fatty acids, especially arachidonic and docosahexaenoic acids; it consumes oxygen at a high rate; and it contains high concentrations of redox-active transition metals but a relatively low concentration of antioxidants. In ALS, at the cellular level, genetic factors, excitotoxicity, apoptosis, inflammation, mitochondrial dysfunction, protein aggregates, and oxidative stress are among the primary hypotheses put forth to explain motor neuron degeneration. Among these factors, oxidative stress appears intimately linked to a series of cellular events in motor neurons that contribute to neuronal degeneration and death [63, 64]. Individuals that had the homozygous genotypes RR and heterozygous GQ had a lower survival rate when compared to the homozygote genotype QQ. Moreover, the allele R was associated with bulbar onset. Some authors have mentioned the increased chance of ALS following the exposure to pesticides and the susceptibility increases two folds with single nucleotide polymorphisms (SNPs) like G-832A, G-162A and C-108 T [65, 66].


4. Conclusions

Distributed watershed models are known as the very powerful tools both for scenario development and for simulating the effects of watershed dynamics management on soil and water resources. This study was aimed to improve understanding of (a) the effects of bypass flows due to irrigation on the calibration of the SWAT model, (b) irrigation return flow and/or drainage generating processes, and (c) N leaching dynamics with simulation of agricultural land management (fertilization, irrigation, and plant species) under the Mediterranean climate conditions. To this aim, the research was conducted in an irrigated agricultural catchment of Akarsu irrigation district. Visual examination of data used in modeling has indicated that drainage flows and nitrogen-leaching processes are not governed by the natural processes in the catchment but mostly by anthropogenic activities.

Model calibration and validation were carried out to determine the most sensitive and appropriate parameter values for the drainage flows generated by the agricultural catchment. Although daily flow data were used in modeling, quantitative model performance evaluation statistics (R2 , NSE, and PBIAS) revealed clearly that the calibrated SWAT model produced rather satisfactory simulation results at the catchment outlet in wet, average, and dry years. In the irrigated catchment, irrigation water losses directly from irrigation channels to drainage ditches, i.e., bypass flows, has direct influence on calibrating hydrologic part of the SWAT model. In this case, the SWAT model findings helped us to highlight that almost 40% of diverted irrigation water has been recklessly squandered in the irrigation scheme. It is almost impossible to quantify bypass flow magnitudes in such irrigation system without using any modeling tools.

Furthermore, modeling exercises showed that the SWAT model run results were sensitive on crop rotations due to the fact that runoff by precipitation and irrigation applications are affected by the land use and land cover types. Contrary to the expectations, daily nitrate modeling results were not able to yield rather satisfactory model performance statistics, indicating that simulated daily nitrogen loads data in drainage were not sufficiently matched with the measured ones. Visual evaluation of measured and simulated nitrogen graphs showed implicit signals that measured nitrogen data might involve some inherent uncertainties and irregularities at the catchment level. Based on the findings, as highlighted in the literature [59], we concluded that model performance can be improved to some extent by increasing the time step from daily to monthly or yearly level for the nitrogen data with involves inherent uncertainties. These uncertainties should be considered when calibrating, validating, and evaluating watershed models because of differences in inherent uncertainty between measured flow, sediment, and nutrient data.

Improved fertilization practices are not only necessary for farmer's economy but also crucial for preserving soil and water resources. In recent years, especial soil analysis in the study area became a very useful tool for fertilizer subsidizes and expert recommendations. However, recommendations can not only be related to and designed by the soil analysis, it should be comprehensively evaluated in a broader environment. At this stage, a suitable model performance enables modeling more sensitive management practices like the fertilizer rates.


*Article* **Mathematical Considerations for Unmanned Aerial Vehicle Navigation in the Magnetic Field of Two Parallel Transmission Lines**

**Dean Martinovi´c, Stjepan Bogdan \* and Zdenko Kovaˇci´c**

Faculty of Electrical Engineering and Computing, University of Zagreb, 10000 Zagreb, Croatia<EMAIL_ADDRESS>(D.M<EMAIL_ADDRESS>(Z.K.)

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** This publication deals with the navigation of unmanned aerial vehicles (UAVs) moving in the magnetic field of two long, straight, parallel conductors, which is of high interest for several new technical applications. How the position and orientation of the UAV can be calculated using a minimal number of only three three-axis magnetometers are discussed. It is shown that the angles can be determined without the knowledge of the conductor currents and the magnetic field equations, but only by combining the sensor measurements with the rotation matrix and exploiting a characteristic property of the magnetic field. Furthermore, different strategies were investigated to determine the respective sensor positions. An analytical solution was derived from the nonlinear magnetic field equations, which promises a low computational time. It is shown that for a given sensor, several solutions exist, from which the correct one has to be selected. Therefore, a specific detection method is introduced. Once the solution is known, the UAV location can be determined. Finally, the overall algorithm was tested by simulations far from and near the conductors with superimposed typical magnetic noise.

**Keywords:** magnetic field navigation; parallel conductors; transmission lines; unmanned aerial vehicles; aerial manipulation

#### **1. Introduction**

Magnetic field-based navigation is an old discipline that had its beginnings in the 20th Century [1]. Magnetic fields have several advantages over other physical quantities. For example, unlike electromagnetic waves, they are not subject to multipath effects or fading because they do not reflect off surfaces. This allows for much higher accuracy. In addition, they are not attenuated when passing through magnetically neutral materials, making magnetic field-based navigation ideal for use between obstacles, at night, or in difficult weather conditions such as snow or fog. A large number of publications exist dealing with techniques for object navigation and localization. Depending on the application, they all use one or more distributed magnetic field sources that emit an alternating magnetic signal. This is usually sampled and measured by several sensors, which may be from different technologies. By combining these measurement data with the mathematical model of the source used, the position and orientation of the observed object are then determined. This can be divided into four different types of sources.

The first type is a simple rectangular or circular coil that typically has a large number of turns. It represents the basic element for the construction of more complex coil types and is often referred to as a dipole, since above a minimum distance from the winding, its magnetic field distribution can be described by the equations of a magnetic dipole. In various research projects [2–6], such as MILPS (magnetic indoor local positioning system), multiple coils mounted on walls are used to locate cell phones and robots in buildings or patients in hospitals. Other use cases include the navigation of medical instruments

**Citation:** Martinovi´c, D.; Bogdan, S.; Kovaˇci´c, Z. Mathematical Considerations for Unmanned Aerial Vehicle Navigation in the Magnetic Field of Two Parallel Transmission Lines. *Appl. Sci.* **2021**, *11*, 3323. https://doi.org/10.3390/app11083323

Academic Editor: Alejandro Suarez

Received: 9 March 2021 Accepted: 2 April 2021 Published: 7 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

inside the human body [7] or the localization of an implantable transponder for intensitymodulated radiation therapy [8]. In addition, several projects have investigated how these coils can be used to locate trapped persons as quickly as possible [9–12]. These could be buried miners or avalanche victims, who can be found by the emitted signal from a small coil worn close to the body. Another interesting topic is the detection of buried landmines in old war zones, also known as humanitarian demining [13,14]. This involves the use of hand-held coils such as metal detectors. Finally, a new application has emerged with the advent of the electric vehicle age.

In the context of inductive charging, the receiver coil mounted on the underbody of the car must be positioned relatively precisely above the transmitter coil in the parking lot, as shown in Figure 1. Since this is very difficult for the driver to achieve, the MPPS (magnetic pulse positioning system) was developed in [15]. Here, the transmitting coil in the ground emits a low-frequency pulsed magnetic signal that allows precise localization by the vehicle in the sub-centimeter range.

**Figure 1.** Navigation using a one-axis coil [15].

The next two coil types consist of two or three base coils wound orthogonally to each other. An example of a three-axis magnetic field source is shown in Figure 2.

**Figure 2.** Principle of a three-axis coil.

For two-axis sources, the coils carry phase quadrature currents that produce a rotating magnetic dipole [1,16,17]. Together with phase detection algorithms and amplitude evaluation, the object position and also the orientation can be determined. For three-axis sources, all source axes are successively excited with currents of the same frequency and phase. In this way, each axis generates an equivalent dipole of arbitrary orientation [18–20]. The three sensor responses are linearly independent and provide sufficient information to determine both the sensor position and its orientation.

The last type of source to be mentioned is a long wire with an alternating signal. In [21–23], its magnetic field was used to allow an underwater vehicle to locate and track a buried underwater cable.

Until now, however, no attention has been paid to magnetic sources consisting of two long parallel conductors. They enable many new useful engineering applications, especially when it comes to the realization of unmanned aerial vehicles (UAVs). For example, UAVs could autonomously perform inventory in warehouses. For this purpose, parallel guides can be attached to the floor or ceiling for navigation between shelves. In addition, UAVs could provide an inspection of rail systems efficiently when a current flows through them.

Another interesting topic is the autonomous performance of inspection and assembly work around power line systems. This is exactly what was investigated in the Aerial-Core (https://aerial-core.eu/, accessed on 06.04.2021) project. One of the goals of the project was to develop a technology that will enable a UAV to perform (Figure 3) the placement of distance separators between power lines.

**Figure 3.** Distance separator for transmission lines.

This is usually performed by workers who are exposed to very high electric and magnetic fields. Moreover, such a task requires a highly skilled helicopter pilot to keep the helicopter close to the power lines while a worker performs the installation. Therefore, the main objective of this work is to derive the equations for determining the position and orientation of a UAV moving through the magnetic field of two parallel conductors, thus enabling precise navigation of a UAV in close proximity to the power lines. We assumed that a UAV is equipped with visual and other sensors that are used both for navigating farther away from the lines and for approaching the lines.

#### **2. Orientation and Localization Equations**

#### *2.1. System Overview*

Figure 4 shows two parallel lines, which together form one phase of a three-phase power transmission system and a UAV flying in close proximity. The transmission lines lie in the *xy*-plane of the global coordinate frame {*G*} positioned exactly between the lines. For the purpose of sensing the magnetic field of the transmission lines, there are magnetometers mounted on the UAV, which are marked with yellow points.

**Figure 4.** System overview. Global frame and UAV frame.

The transmission lines of a single phase carry sinusoidal currents of fixed frequency and phase, which can be used to our advantage. They generate an alternating magnetic field, which allows static and low-frequency magnetic interference fields such as the Earth's field to be canceled out. This is realized by a differential measurement of the amplitude of the magnetic signal. The calculation can be implemented in the time domain and is very efficient [24]. In addition, the system can be treated as a static field problem with constant currents. For simplicity, they are chosen as *I* in the following. The magnetic field thus generated can be calculated using the law of Biot–Savart for thin conductors. In general, transmission lines are not thin. However, as the distance from the lines increases, the magnetic flux density distribution approaches that of thin conductors relatively quickly:

$$\vec{B}(\vec{r}) = -\frac{\mu\_0}{4\pi} I \int\_{-\infty}^{\infty} \frac{(\vec{r} - \vec{r}\,') \times d\vec{s}\,'}{|\vec{r} - \vec{r}\,'|^3} \,. \tag{1}$$

The integration is executed only along the cables according to the arrangement in Figure 4, i.e., *ds* = (*dx* , 0, 0)*T*. Thus, from the cross-product in Equation (1), it immediately follows that:

$$B\_{\mathbf{x}} = \mathbf{0}, \forall (\mathbf{x}, \mathbf{y}, \mathbf{z})^{T}, \tag{2}$$

which means that there is no field component in the direction parallel to the cables. Solving the integral for one cable and using the superposition principle, the remaining field components *By* and *Bz* can be determined. The total magnetic field equation with the transmission lines at ±*y*0 is then:

$$\vec{B}(y, z) = \mathcal{C} \begin{pmatrix} 0\\ \frac{-z}{z^2 + (y + y\_0)^2} + \frac{-z}{z^2 + (y - y\_0)^2} \\ \frac{y + y\_0}{z^2 + (y + y\_0)^2} + \frac{y - y\_0}{z^2 + (y - y\_0)^2} \end{pmatrix}, \quad \text{with} \quad \mathbb{C} = \frac{I\mu\_0}{2\pi}, \tag{3}$$

where *C* is a prefactor depending on the transmission line current.

Figure 5 shows the signal amplitude |- *B*(*y*, *z*)| around the transmission lines for *I* = 700 A and *y*0 = 0.25 m. Figure 6 shows the locations of equal amplitude in the *yz*-plane expressed by isolines. The magnetic field contributions of the other two phases are neglected because the cables are typically at least 2 m apart.

**Figure 5.** Signal power around the transmission lines.

**Figure 6.** Isolines of the magnetic field describing equal power.

The magnetic field-based navigation of a UAV requires that the UAV can sense the magnetic field in all spatial directions. It uses this information to determine its relative position and orientation with respect to the fixed position and orientation of the global coordinate frame {*G*}. The flight capabilities allow the UAV to achieve an arbitrary orientation by rotating the UAV's body about the *z*- (blue), *y*- (green), and *x*-axes (red) of the local drone frame *Lu* for the respective yaw, pitch, and roll angles *α*, *β*, and *γ*. The origin position of the local frame *Lu* in {*G*}, denoted by *pD*, also represents the position of the UAV to be determined. The sensors have predefined positions *pi*0 in *Lu*, where the index *i* refers to a particular sensor. Moreover, each magnetometer has its own local coordinate frame *Li*. Then, their respective positions *pi* = (*xi*, *yi*, *zi*)*T* in {*G*} can be expressed using the rotation matrix *RLu G* (*α*, *β*, *γ*) = *Rx*(*γ*) *Ry*(*β*) *Rz*(*α*):

$$
\vec{p}\_{\vec{l}} = \mathcal{R}\_{G}^{L\_{\text{av}}}(\mathfrak{a}, \mathfrak{z}, \mathfrak{z}) \,\vec{p}\_{\vec{n}0} + \vec{p}\_{D\prime} \,\tag{4}
$$

$$R\_{\mathbb{Z}}(\mathfrak{a}) = \begin{pmatrix} c\_{\mathfrak{a}} & s\_{\mathfrak{a}} & 0 \\ -s\_{\mathfrak{a}} & c\_{\mathfrak{a}} & 0 \\ 0 & 0 & 1 \end{pmatrix}, \\ R\_{\mathbb{Y}}(\mathfrak{f}) = \begin{pmatrix} c\_{\mathfrak{f}} & 0 & -s\_{\mathfrak{f}} \\ 0 & 1 & 0 \\ s\_{\mathfrak{f}} & 0 & c\_{\mathfrak{f}} \end{pmatrix}, \\ R\_{\mathbb{X}}(\gamma) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & c\_{\gamma} & s\_{\gamma} \\ 0 & -s\_{\gamma} & c\_{\gamma} \end{pmatrix}.$$

where *sx* and *cx* are abbreviations for *sin x* and *cos x* needed in the following. From Equation (3), it can be seen that the position of the UAV along the cables (*x*-direction) cannot be determined because the magnetic field does not provide sufficient information. Thus, only the projections of the magnetometers and the UAV frame origin in the *yz*-plane can be computed. The representation of the magnetic field *vi* measured by the *i*-th sensor at position *pi* expressed in {*G*} then becomes:

$$\vec{B}(\vec{p}\_i) = \mathcal{R}\_G^{L\_u}(\alpha, \mathcal{G}, \gamma) \; R\_{L\_u}^{L\_i} \; \vec{v}\_{i\prime} \tag{5}$$

where *RLi Lu* is the user-defined orientation of *Li* in the UAV body frame *Lu* and *RLu G* is the orientation of *Lu* in the global frame {*G*}. If *Li* is aligned with {*G*}, then it is *RLu G RLi Lu* = *I*, where *I* is the identity matrix. In this case, the sensor measures exactly the magnetic field given by Equation (3), i.e., - *B*(*pi*) = *vi*. The system of three equations in (5) establishes a relationship between the sensor measurement, sensor position, and UAV orientation. It is nonlinear and has five unknowns (*yi*, *zi*, *α*, *β*, *γ*)*T*, which are determined in the following sections.

#### *2.2. Orientation Determination*

In this section, the UAV angles *α*, *β*, and *γ* are derived. As will be shown below, the calculation is possible without knowledge of the field equation and the transmission

line current. Three magnetometers were considered. For the following analysis, their corresponding equations were set up in matrix notation, where 03×3 is a 3 × 3 zero matrix:

$$
\begin{pmatrix}
\vec{B}(\vec{p}\_1) \\
\vec{B}(\vec{p}\_2) \\
\vec{B}(\vec{p}\_3)
\end{pmatrix} = \begin{pmatrix}
R\_G^{L\_u} R\_{L\_u}^{L\_1} & 0\_{3\times 3} & 0\_{3\times 3} \\
0\_{3\times 3} & R\_G^{L\_u} R\_{L\_u}^{L\_2} & 0\_{3\times 3} \\
0\_{3\times 3} & 0\_{3\times 3} & R\_G^{L\_u} R\_{L\_u}^{L\_3}
\end{pmatrix} \begin{pmatrix}
\vec{v}\_1 \\
\vec{v}\_2 \\
\vec{v}\_3
\end{pmatrix}.\tag{6}
$$

It is obvious that this system of equations is solvable since the rank of the matrix is nine, which is equal to the number of unknowns. In general, in addition to the three angles, each row provides two unknowns *yi* and *zi*. Now, the UAV's yaw and pitch angles can be calculated directly from the two sensor measurements *vl* and *vn* and the zero field property (2). With this and the Cartesian unit vector*ex*, it follows from Equation (6) that:

$$0 = (R\_G^{L\_u} R\_{L\_u}^{L\_i} \vec{v}\_i) \cdot \vec{e}\_{x\_{\prime}} \qquad i = l, n. \tag{7}$$

Obviously, this equation does not depend on the sensor positions or the choice of their local frames. Thus, they can be located anywhere on the UAV and have any local frame. However, to avoid additional transformations or computational effort, the sensor frames are aligned with the UAV body frame. Then, the equation system (7) becomes:

$$\begin{array}{ccccc} 0 = & c\_{\beta} \ c\_{a} \ v\_{lx} & + c\_{\beta} \ s\_{a} \ v\_{ly} & -s\_{\beta} \ v\_{lz} \\ 0 = & c\_{\beta} \ c\_{a} \ v\_{nx} & + c\_{\beta} \ s\_{a} \ v\_{ny} & -s\_{\beta} \ v\_{nz} \end{array} \tag{8}$$

After division by *cβ* and subsequent subtraction, *tan β* can be eliminated. Solving the remaining equation for *α* leads to:

$$\alpha = \operatorname{atan} \left( \frac{\upsilon\_{lz}\upsilon\_{nx} - \upsilon\_{lx}\upsilon\_{nx}}{\upsilon\_{ly}\upsilon\_{nz} - \upsilon\_{lz}\upsilon\_{ny}} \right) = \operatorname{atan} \left( \frac{(\vec{\upsilon}\_{l} \times \vec{\upsilon}\_{n}) \cdot \vec{\varepsilon}\_{y}}{(\vec{\upsilon}\_{l} \times \vec{\upsilon}\_{n}) \cdot \vec{\varepsilon}\_{x}} \right) \tag{9}$$

and further with *A* = *tan α* to:

$$\beta = \operatorname{atan} \left( \frac{v\_{nx} + v\_{ny}A}{v\_{nz}\sqrt{A^2 + 1}} \right) = \operatorname{atan} \left( \frac{v\_{lx} + v\_{ly}A}{v\_{lx}\sqrt{A^2 + 1}} \right). \tag{10}$$

The roll angle *γ* can also be computed by exploiting the zero field property (2). However, the rotation matrix *RLu G* does not in general contain *γ* in its *x*-row. Thus, a different sensor *k* must have a different local frame *Lk* in order for the zero to appear in the *y*- or *z*-row containing all three angles. This can be achieved by, for example, placing the sensor so that it is rotated relative to the UAV body frame *Lu* by an angle *δ* about the *y*-axis, i.e., it is:

$$\mathcal{R}\_{L\_{\mu}}^{L\_{k}} = \mathcal{R}\_{\mathcal{Y}}(\boldsymbol{\delta}).\tag{11}$$

Then, from Equation (5), it follows that after applying the coordinate transformation *A* = *RLu G RT y* (*δ*)*RG Lu* , the vector - *B*(*pk*) is expressed in *Lu*:

$$\begin{split} \vec{\vec{B}}(\vec{p}\_k) &= A \, \vec{B}(\vec{p}\_k) = R\_G^{L\_w} R\_\mathcal{Y}^T(\delta) R\_{L\_w}^G R\_G^{L\_w} R\_\mathcal{Y}(\delta) \, \vec{v}\_{k\prime} \\ &\quad \text{for } \delta = \frac{\pi}{2} \Rightarrow \begin{pmatrix} \mathcal{B}\_x(\vec{p}\_k) \\ \hat{\mathcal{B}}\_\mathcal{Y}(\vec{p}\_k) \\ 0 \end{pmatrix} = R\_G^{L\_w} \vec{v}\_k. \end{split} \tag{12}$$

From Equation (12), it follows that the rotation of - *B*(*pk*) for a given angle *δ* = *π* 2 around the *y*-axis of *Lu* has forced the *z*-component of the resulting vector - *B*ˆ(*pk*) to zero:

$$\mathcal{B}\_z(\vec{p}\_k) = 0 = (R\_G^{L\_u} \vec{v}\_k) \cdot \vec{e}\_{z\prime} \tag{13}$$

where*ez* is the Cartesian unit vector. Then, in a first step, Equation (13) can be written as:

$$0 = (s\_{\gamma}s\_{a} + c\_{\gamma}c\_{a}s\_{\beta})\upsilon\_{kx} + (c\_{\gamma}s\_{\beta}s\_{a} - c\_{a}s\_{\gamma})\upsilon\_{ky} + c\_{\gamma}c\_{\beta}\upsilon\_{kz}.$$

Next, the division by *cγ* forms the tangent function *tγ*, which can be isolated and finally solved for *γ*:

$$\begin{split} 0 &= (t\_{\gamma}s\_{\mathfrak{a}} + c\_{\mathfrak{a}}s\_{\beta})v\_{kx} + (s\_{\beta}s\_{\mathfrak{a}} - c\_{\mathfrak{a}}t\_{\gamma})v\_{ky} + c\_{\beta}v\_{kx} \\ &\Rightarrow \gamma = \operatorname{atan}\left(\frac{v\_{ky}s\_{\beta}s\_{\mathfrak{a}} + v\_{kx}s\_{\beta}c\_{\mathfrak{a}} + v\_{kx}c\_{\beta}}{v\_{ky}c\_{\mathfrak{a}} - v\_{kx}s\_{\mathfrak{a}}}\right) . \end{split} \tag{14}$$

The sensor can be placed anywhere on the UAV. However, the result only applies to the local frame definition in (11). Theoretically, arbitrary local frame definitions can be chosen, but in this case, there is no zero in any coordinate, and the whole system of Equation (12) must be solved since the sensor position *pk* is unknown. However, the equations are complicated, and it is questionable whether an analytical solution exists. The presented formulas are valid if the UAV angles are within the definition range *D*, respectively *α*, *β*, *γ* ∈ *D* = [−*π* 2 , *π* 2 ]. The limits of the equations presented are discussed in Section 3.

#### *2.3. Position Determination*

In the previous section, it was shown that the UAV orientation and position computation can be decoupled into two independent problems. In doing so, it was shown how the UAV orientation can be computed using the zero field property and without knowledge of the transmission line currents or the magnetic field equation. This section now deals with the computation of the respective sensor positions *pi*, which are needed to finally determine the UAV position *pD* via Equation (4).

With the decoupling presented, the complexity of the equation system (6) can be greatly reduced, since by knowing the angles, the right side can be computed and becomes a simple real-numbered vector. With *RLu G RLi Lu vi* = *Vi*, the system becomes:

$$
\begin{pmatrix}
\vec{B}(\vec{p}\_1) \\
\vec{B}(\vec{p}\_2) \\
\vec{B}(\vec{p}\_3)
\end{pmatrix} = \begin{pmatrix}
\vec{\mathbb{V}}\_1 \\
\vec{\mathbb{V}}\_2 \\
\vec{\mathbb{V}}\_3
\end{pmatrix}.
\tag{15}
$$

Obviously, the computation of the respective sensor positions is also decoupled. Thus, to obtain *pD*, it is sufficient to solve one of these rows in (15). Then, for a sensor *i*, the corresponding row is:

$$
\vec{B}(\vec{p}\_i) = \begin{pmatrix} B\_y(y\_{i\prime}z\_i) \\ B\_z(y\_{i\prime}z\_i) \end{pmatrix} = \begin{pmatrix} V\_{\rm iy} \\ V\_{\rm iz} \end{pmatrix} = \vec{V}\_{\rm i\prime} \tag{16}
$$

which represents a system of two equations. The components *By* and *Bz* describe the magnetic field in two dimensions in {*G*} and were obtained by calculating the Biot–Savart integral (1). This means, on the one hand, that the two equations in (16) are linearly independent and theoretically allow the determination of the unknowns *yi* and *zi*. On the other hand, the equations are nonlinear, and the question arises whether it is possible to solve (16) analytically and, if so, how. In this context, there are several approaches, all of which were investigated in the context of this work. Since the following considerations are equally valid for all sensors, the subscript *i* is omitted for simplicity.

The first row in (16) can be solved for one coordinate *y* or *z*, which is then inserted into the second row to compute the other coordinate; and the same in reverse. However, solving *By* = *Vy* for *z* or *Bz* = *Vz* for *y* in the general form leads to extremely large solutions *z*(*y*, *Vy*) and *y*(*z*, *Vz*), which contain root terms with eighth-order polynomials as arguments. The

other two cases, i.e., solving *By* = *Vy* for *y* or *Bz* = *Vz* for *z*, lead to smaller solutions. For example, from the latter, the four solutions for *z* follow:

$$z(y, V\_z) = \pm \sqrt{\pm D(y^2, y) + a(y^2, y)},$$

$$\begin{array}{ll}\text{with} & D = \frac{1}{V\_z}\sqrt{4V\_z^2 y^2 y\_0^2 - 4CV\_z y y\_0^2 + C^2 y^2} \\ & a = -y\_0^2 - y^2 + \frac{Cy}{V\_z} \end{array} \tag{17}$$

The result *y*(*z*, *Vy*) has the same structure, but is not used, as the easiest final equation can be found by inserting (17) into *B*2 *y*(*y*, *z*) = *V*2 *y* . With *D* from (17), it becomes:

$$\frac{V\_z(\text{Cy} \pm D)^2 (-V\_z y\_0^2 - V\_z y^2 + \text{Cy} \pm D)}{y^2 (-2V\_z y\_0^2 + \text{Cy} \pm D)^2} = V\_y^2. \tag{18}$$

The equation can in principle be solved for *y*. However, no suitable substitution could be found, so that the root terms disappear, and at the same time, a polynomial of a smaller order than five is obtained. For example, after the substitution with *E* = *Cy* ± *D*, unwanted *y*2 expressions still remain in the equation. The square roots can also be eliminated by subsequent isolation and squaring. However, this leads to a high-order polynomial. In summary, it is difficult to impossible to derive an analytical solution only from the equation system (16). In this context, the idea is to consider the signal power as an additional equation. As will be shown in the next steps, its combination with (16) leads to a compact analytical solution, which is derived below: At a given position, the signal power can be calculated by taking the dot product of the corresponding sensor response *V* with itself. This must be equal to the dot product of - *B*(*y*, *z*) with itself, leading to the following equation:

$$|\vec{V} \cdot \vec{V} = |\vec{V}|^2 = \vec{B}(y, z) \cdot \vec{B}(y, z) \Rightarrow \quad |\vec{V}|^2 = \frac{4\mathcal{C}^2(z^2 + y^2)}{z^4 + y^4 + (2y\_0^2 + 2y^2)z^2 + y\_0^4 - 2y^2y\_0^2}.\tag{19}$$

The most simple equation is found by solving the formula for *z*2 and inserting it into *Bz*(*y*, *z*) = *Vz*. For this, the two solutions for *z*2 are obtained with *P* = |*V*-|:

$$z^2(y, P) = \pm 2\sqrt{y^2 y\_0^2 - \frac{C^2 y\_0^2}{P^2} + \frac{C^4}{P^4}} - y\_0^2 - y^2 + \frac{2C^2}{P^2}.\tag{20}$$

Next, with the positive solution of (20), the equation *Bz*(*y*, *z*2(*y*, *P*)) = *Vz* can be written as:

$$\frac{-P^4 y g\_0^2 + \hat{D} P^2 y + \mathcal{C}^2 P^2 y}{-\mathcal{C} P^2 y\_0^2 + 2\mathcal{C}\hat{D} + 2\mathcal{C}^3} = V\_{z\_{\prime}} \quad \text{with} \quad \hat{D} = \sqrt{(P^4 y^2 - \mathcal{C}^2 P^2) y\_0^2 + \mathcal{C}^4}.$$

Further, *D*ˆ is isolated and squared, which leads to:

$$D^2 = \frac{((P^4y - CP^2V\_z)y\_0^2 - C^2P^2y + 2C^3V\_z)^2}{(P^2y - 2CV\_z)^2}.$$

This eliminates the square root of *D*ˆ . After re-substituting and transforming the equation, the position calculation is reduced to a root-finding problem of a fourth-order polynomial,

$$\begin{array}{cccc}\sum\_{n=0}^{4}a\_{n}y^{n} = 0, & a\_{4} = P^{8}y\_{0\prime}^{2} & a\_{3} = -4\mathbb{C}P^{6}V\_{z}y\_{0\prime}^{2} \\ a\_{2} = -P^{8}y\_{0}^{4} + \mathbb{C}^{2}P^{6}y\_{0}^{2} + 4\mathbb{C}^{2}P^{4}V\_{z}^{2}y\_{0\prime}^{2} \\ a\_{1} = 2\mathbb{C}P^{6}V\_{z}y\_{0}^{4} - 2\mathbb{C}^{3}P^{4}V\_{z}y\_{0\prime}^{2} & a\_{0} = -\mathbb{C}^{2}P^{4}V\_{z}^{2}y\_{0\prime}^{4}.\end{array} \tag{21}$$

The general analytic root expressions of such a polynomial are extremely large. Therefore, it is important to make Equation (3) as simple as possible by placing the global coordinate frame exactly halfway between the transmission lines. In addition, the lines

should be placed in the *xy*-plane parallel to the *x*-axis. Then, with the coefficients *an* given in (21), the four large expressions for *y* are greatly reduced and become:

$$y = \pm \frac{\sqrt{P^4 y\_0^2 + 2C^2 V\_z \sqrt{V\_z^2 - P^2} + 2C^2 V\_z^2 - C^2 P^2}}{2P^2} \pm \frac{\sqrt{P^4 y\_0^2 - 2C^2 V\_z \sqrt{V\_z^2 - P^2} + 2C^2 V\_z^2 - C^2 P^2}}{2P^2} + \frac{C V\_z}{P^2} \tag{22}$$

where *S* represents the two square root addends. The same result is obtained for the negative solution of Equation (20). Now, in the next step, four different cases can be distinguished, all of which lead to a certain class of solutions.

(I) *Vy*, *Vz* = 0:

**Theorem 1.** *Two of the y-solutions in* (22) *can immediately be excluded, because they are complex.*

**Proof.** *Vz* ≤ |*V*- | = *P*, ∀*V*- ∈ R3 ⇒ &*V*2 *z* − *P*2 ∈ C. With this, the above square root addends represented by *S* can be written as:

$$S = \pm \sqrt{a + jb} \pm \sqrt{a - jb},$$

$$a = P^4 y\_0^2 + (2V\_z^2 - P^2)C^2, \; b = 2C^2 V\_z \sqrt{P^2 - V\_z^2}.\tag{23}$$

Next, the parameters *r* and *ϕ* are introduced as:

$$
\sigma = \sqrt{a^2 + b^2}, \quad \varrho = \tan^{-1}\left(\frac{b}{a}\right) + \pi. \tag{24}
$$

Now, if the addends have a different sign *S*, this results in:

$$\begin{split} S &= \pm \sqrt{a+j}\,b \mp \sqrt{a-j}\,b = \pm \sqrt{r}\,e^{j\frac{\theta}{2}} \mp \sqrt{r}\,e^{-j\frac{\theta}{2}}, \\ &= \pm 2j\sqrt{r}\,\frac{1}{2\overline{\gamma}}(e^{j\frac{\theta}{2}} - e^{-j\frac{\theta}{2}}) = \pm 2j\sqrt{r}\sin\frac{\theta}{2}, \end{split}$$

whereby the *e*-expressions are obtained using Euler's formula.

Applying this procedure for the case of the same sign, the real-numbered and more compact expressions for *y* can be derived. Therefore, in a first step, *S* is determined:

$$S = \pm \sqrt{a+j'b} \pm \sqrt{a-j'b} = \pm \sqrt{r} \, e^{j\frac{\theta}{2}} \pm \sqrt{r} \, e^{-j\frac{\theta}{2}}$$

$$= \pm 2\sqrt{r} \, \frac{1}{2} (e^{j\frac{\theta}{2}} + e^{-j\frac{\theta}{2}}) = \pm 2\sqrt{r} \cos \frac{\theta}{2}.$$

The square root addends in Equation (22) represented by *S* now can be replaced by this result, leading to:

$$y\_{1,2} = \frac{S}{2P^2} + \frac{CV\_\tau}{P^2} = \frac{1}{P^2} \left( \pm \sqrt{r} \cos \frac{\varphi}{2} + CV\_z \right),\tag{25}$$

with *r* and *ϕ* according to (24).

**Theorem 2.** *There exist two solutions* |*y*1| = |*y*2|*. One each is in y* < 0 *and y* > 0*, respectively.*

**Proof.** The hypothesis is true if the absolute value of the square root addends in Equation (22) is always greater than *C*|*Vz*|*P*−2.

$$\begin{split} \frac{|\pm(\sqrt{a+jb}+\sqrt{a-jb})|}{2P^2} > \frac{\mathbb{C}|V\_z|}{P^2} &\implies |(\sqrt{a+jb}+\sqrt{a-jb})|^2 > (2\mathbb{C}|V\_z|)^2 \\ \Rightarrow q\_l &= \sqrt{a+jb}\sqrt{a-jb} > 2\mathbb{C}^2V\_z^2 - a = q\_r. \end{split}$$

Further, with *a*, *b* from (23), it follows for the left and right side:

$$q\_I = -P^8 y\_0^4 + (4\mathcal{C}^2 P^4 V\_z^2 - 2\mathcal{C}^2 P^6) y\_0^2 + \mathcal{C}^4 P^4,$$

$$q\_I = -P^8 y\_0^4 - 2\mathcal{C}^2 P^6 y\_0^2 + \mathcal{C}^4 P^4,$$

$$\Rightarrow q\_I - q\_I = 4\mathcal{C}^2 P^4 V\_z^2 y\_0^2 > 0.$$

(II) *Vy* = 0 ⇒ *P* = |*Vz*|: In this case, the excluded complex solutions must also be considered, since they become real. Both give the same result *y*3. The two real solutions *y*1,2 can again be obtained from (25) or from the general Equation (22).

$$y\_{1,2} = \frac{1}{V\_z^2} (\pm \sqrt{V\_z^4 y\_0^2 + \mathcal{C}^2 V\_z^2} + \mathcal{C}V\_z), \quad y\_3 = \frac{\mathcal{C}}{\nabla\_z} \tag{26}$$

(III) *Vz* = 0 ⇒ *P* = |*Vy*|: The solutions are obtained analogous to Case (II). |*y*1| = |*y*2|.

$$y\_{1,2} = \pm \frac{\sqrt{V\_y^4 y\_0^2 - \mathcal{C}^2 V\_y^2}}{V\_y^2}, \quad y\_3 = 0 \tag{27}$$

(IV) *Vy* = *Vz* = 0: If this case occurs, the sensor lies exactly in the origin of the global frame.

Once the corresponding *y* is determined, in a final step, *z* can be calculated from Equation (20):

$$z = \pm \sqrt{\pm 2\sqrt{(y^2 - \frac{C^2}{P^2})y\_0^2 + \frac{C^4}{P^4}} - y\_0^2 - y^2 + \frac{2C^2}{P^2}},\tag{28}$$

offering four possible *z*-solutions, all of which can be real-numbered. There are also situations where, for a given *y*, a subset of the *z*-solutions in (28) can become complex. Therefore, their use must be treated with caution, especially when the formulas are implemented on a microcontroller of a real UAV. A safer way is to check the arguments of the square roots. The argument of the outer square root can be written as:

$$A\_o = \pm g + h, \; h = -y\_0^2 - y^2 + \frac{2C^2}{P^2}, \tag{29}$$

where *g* is the inner square root expression. For all solutions *y*1,2,3 obtained, *g* ≥ 0 holds, since they were obtained using the *z*2-expression (20). Thus, the argument of the inner square root does not need to be checked. In contrast, *Ao* can become negative in general and should be checked to satisfy that the condition *Ao* > 0 is fulfilled before applying the square root.

#### *2.4. Selecting the Right Solution*

Equation (28) yields four *z*-solutions for each *y*-solution of Cases (I)–(IV). Thus, depending on the case, there are up to twelve possible positions in the *yz*-plane for a given sensor. For navigation, it is necessary to know which of these points is the correct one. Finding it is the task of this subsection.

Now, in the first approach, the number of possible positions can be narrowed down by excluding Cases (II)–(IV). The reason for this is that they can only occur in very few locations and for a very short time during navigation, which extremely reduces their probability of occurrence. For example, the solutions of Case (II) all lie on the *z*-axis. In Case (III), one solution lies on the *y*-axis and the other two lie on a circle with a diameter equal to the distance between the transmission lines, i.e., 2*y*0. This can be shown by computing

the radius ' *y*2 1,2 + |*z*2(*y*1,2)|, where *y*1,2 comes from Equation (27) and *z*2 comes from (20). The case with the lowest probability of occurrence is Case (IV), where the sensor is located directly at the origin. The locations of Cases (II)–(IV) are summarized in Figure 7 and shown as black and red lines. Since they are not relevant for navigation and only increase the implementation effort, they can be ignored. Thus, since only Case (I) is considered, eight possible positions for a single sensor remain, from which the correct one must still be selected.

In the next steps, it will be seen that six can be excluded. To do this, some constraints must be met. As explained in Section 2.1, the sinusoidal transmission line current can be considered as an equivalent constant current *I*, i.e., it never changes its direction. Furthermore, the UAV orientation must not exceed the definition range *D* from Section 2.2 with *α*, *β*, *γ* ∈ *D* = [−*π* 2 , *π* 2 ]. This would otherwise lead to ambiguities due to the symmetry of the magnetic field. The details concerning ambiguities are the subject of Section 3. Now, the idea is to investigate which of the eight remaining positions lead to a match with the sensor measurement *V*-, i.e., which positions satisfy Equation (16).

**Theorem 3.** *Given a sensor k that measures Vk =* (*Vky*, *Vkz*)*T and has two y-solutions according to* (25) *in Case (I), for each of the two values, there are four associated z-solutions according to Equation* (28)*. However, for each y, there is only one z-solution that leads to a match with Vk, i.e., there are only two possible locations in the yz-plane where the sensor can lie.*

**Proof.** Let *yk* be one of the two *y*-solutions. Then, it generates the four *z*-solutions, which in simplified notation are *zn* = ±&±*g* + *h*, *g*, *h* = 0. *zn* has two outer and two inner signs. Further, the *Bz*-component in Equation (3) has a *z*2-addend in the nominator as the only *z*-dependency, so it can be written as *Bz*(*yk*, *z*2). Now, if for *zl* = ±&*g* + *h*, it follows that:

$$B\_z(y\_{k'}z\_l^2) = B\_z(y\_{k'} \mathfrak{g} + h) = V\_{kz'}$$

then with *g* + *h* = −*g* + *h*:

$$V\_{kz} = B\_z(y\_{k\prime} \text{g} + h) \not\equiv B\_z(y\_{k\prime} - \text{g} + h),$$

because *Bz*(*yk*, *z*2) is strictly monotonous. This means, that only one of the inner signs of *zn* leads to a match with *Vkz*. This fact is independent of the outer sign. Once the inner sign is known, similarly, the *By*-component in Equation (3) can be analyzed. Therefore, let *z*2 *m* = *g* + *h* lead to a match with *Vkz*. Now, if *By*(*yk*, *zm*) = *Vky*, then *By*(*yk*, −*zm*) is:

$$\mathcal{L}\left(\frac{-(-z\_m)}{z\_m^2 + (y\_k + y\_0)^2} + \frac{-(-z\_m)}{z\_m^2 + (y\_k - y\_0)^2}\right) = -V\_{ky}.\tag{30}$$

As shown in the last step, only one of the outer signs of *zn* can be the correct one, leading to a match with *Vky*.

In summary, there are only two possible locations for a single sensor, so that Equation (16) is satisfied. With the constraints defined above, another fact follows from the proof. It is needed for the consideration of the ambiguities in Section 3.

**Corollary 1.** *The sign of Vy of the sensor response V with respect to the global frame* {*G*} *depends only on the half plane z* > 0 *or z* < 0 *in which the sensor is located.*

Consequently, both solutions are always in the same half plane. By Theorem 2, one lies in *y* < 0 and the other in *y* > 0. An example is given in Figure 7, which shows the solutions as bold dots.

**Figure 7.** Two remaining positions of a sensor.

Since both points must be in areas with the same direction of the magnetic field lines to satisfy Equation (16), one is outside the circle and one is inside. Of the two remaining sensor positions that match *V*- , one must still be excluded. By considering only a single sensor, it is impossible to determine the correct one. However, in the following, multiple sensors are considered instead, which provide much more information through which a whole set of correct sensor positions can be selected.

The idea is as follows. On the one hand, the user-defined positions *pi*0 of *N* sensors form an original polygon in the *yz*-plane of the UAV frame *Lu*. It has a fixed geometric structure *Su* like a fingerprint. On the other hand, the two positions *pi*1,2 computed for each sensor can now be combined into 2*N* polygons in the global frame {*G*}. Finally, the one that matches the structure *Su* is selected. Therefore, before comparison, each polygon must be transformed back to *Lu* using Equation (4). However, this is not possible since the UAV position *pD* must already be known for this. One way around this problem is to choose a different reference point. For example, all polygons are shifted so that their upper left corner appears in the origin of {*G*} before comparison. What is important here is that the original polygon is first transformed to {*G*}. Its new geometric structure *SG* can be computed using the known UAV orientation and *pi*0 via Equation (4) with *pD* =-0. Then, its upper left corner is moved to the origin. The other polygons obtained from the sensor measurements are already representations with respect to {*G*}, so their upper left corner only needs to be shifted to the origin. The shift is allowed since it does not affect the polygon geometry. Finally, the comparison can be done using the least- squares estimation. Here, the deviation of the corresponding polygon nodes is squared and added up. The polygon with the smallest sum best fits *SG* and is selected. From each of its nodes, the UAV position *pD* can then finally be computed using Equation (4). It is recommended to place the sensors in the *yz*-plane of *Lu* so that they span an area as large as possible. This allows for better polygon detection, since during the aerial manipulation, the roll and pitch angle of the UAV will typically have a rather small misalignment with respect to the global frame {*G*}. The polygon detection method proposed here was successfully tested. The implementation details and test results are presented in Section 4.

#### **3. Ambiguities**

#### *3.1. Ambiguities Due to Field Symmetry*

An ambiguity exists if all of the UAV's sensors show the same output in at least two different locations. Under the condition of equivalent static transmission line current, several cases can be identified, which are shown in Figure 8.

**Figure 8.** Ambiguities due to field symmetry.

In the lower left, the real position of the local frame *Li* of the *i*-th sensor is shown. The UAV can now change its flight direction and additionally be above the transmission line in Case (a). The sensor output would be the same. In Case (b), the whole UAV is rotated by angle *π* around the *x*-axis. Here, since the magnetic field is mapped to itself, the sensor will again measure the same. In this case, the UAV is flying on its back.

#### *3.2. Ambiguities Due to Signal Symmetry*

As introduced in Section 2.1, transmission lines carry an alternating current. In addition to improved filtering capabilities, this also introduces a number of ambiguities in navigation, as shown below. The magnitude and sign of the three field components measured by a three-axis magnetometer can be determined in the time domain from the samples. For this purpose, a method was proposed in [24], which is beyond the scope of the present paper. It is taken as given for the following considerations. For example, a sensor with the sample values depicted in Figure 9 would yield the measured value *v* = (1 *mT*, −200 μ*T*, −500 μ*T*)*T*.

**Figure 9.** Sample values of a magnetometer.

However, a half period later, the signals in Figure 9 appear shifted, and the signs of the individual field components change; i.e., the sensor now measures −*v*. Since technically, only signal amplitudes are used, it is as if the equivalent static current *I* changes its direction every half period. A direct consequence of this fact follows from Corollary 1. Under the defined constraints in Section 2.4, it allows detecting whether the UAV is flying under or above the transmission lines, since the sign of the measured *y*-components *Viy* in the global frame depends only on which half plane *z* > 0 or *z* < 0 the sensors are located. However, since the direction of *I* changes cyclically, this is no longer possible without additional detection methods, since *Viy* also changes sign. Such a detection could be realized, for example, by comparing the signal energy of several sensors. This could exploit the fact that the closer the sensor is to the cables, the higher the energy is. However, this is beyond the scope of this paper. The fact that it is impossible to uniquely determine the sensor

position can even be seen directly from Equation (25). It contains the prefactor *C*, which, as shown in (17), depends on the current; thus, it can be *C* > 0 or *C* < 0, leading to different solutions. Furthermore, the whole approach proposed in Theorem 3 to narrow down the solutions no longer works.

#### *3.3. Resolving the Ambiguities*

If no additional electronic means are to be used, the ambiguities can be resolved as suggested in [15], where an electric vehicle approaches the inductive charging coil of a parking lot. In the case of using only two magnetometers, it never knows whether the coil is in front or behind it. The vehicle therefore assumes that the coil is always in front of it or that the driver is always moving ahead of the parking space. In the same way, concrete constraints for the navigation of a UAV must be defined and passed to it as a priori knowledge. These are:

Furthermore, a selection criterion for the measured values must be established. This is because, as mentioned earlier, the sign of the measurement *vi* of the sensor *i* changes with each half period. It must now be decided which is the correct one. The selection criterion can be derived from the above constraints. With *C* > 0 and *z* < 0 from Conditions (1) and (2), it follows from the right-hand rule that the magnetic field always has a positive *y*component, i.e., *By*(*y*, *z*) > 0. This can also be shown directly with Equation (3). Moreover, according to Theorem 3, this circumstance never changes as long as the UAV remains in the half plane *z* < 0 and as long as Constraint (3) is not broken. Using Equation (5) and the Cartesian unit vector*ey* leads to the necessary selection criterion:

$$(R \; \vec{v}\_i) \cdot \vec{e}\_y > 0. \tag{31}$$

If *vi* does not satisfy (31), −*vi* does, and vice versa. Only one sign can be a match. Direct comparison of the signs of *viy* and *By*(*yi*, *zi*) is not allowed, since from *By*(*yi*, *zi*) > 0, *viy* > 0 does not necessarily follow. The last point to be clarified is the orientation of the UAV. Since the rotation matrix in (31) is used, it must be ensured that the angles can be computed independently of the direction of *I*.

**Theorem 4.** *Given two sensors i and k with equally aligned local frames and fixed locations pi* = *pk, as long as their measurement values vi and vk are obtained from the same signal period, the yaw angle α and pitch angle β are independent of the direction of the static current I.*

**Proof.** Referenced to the same time interval containing a whole signal period, the sensors measure *vi*,*vk* =-0. Half a signal period later, when both measurements are calculated from the same new time interval and the sensors have not moved, the sign of the measurements is reversed. Then, the cross-product in Equation (9) for the yaw angle becomes −*vi* × −*vk*. Now, by applying the rule (*ra*) ×*b* =*a* × (*r b*) = *r*(*a* ×*b*), it follows that:

$$-\vec{\upsilon}\_i \times -\vec{\upsilon}\_k = (-1)^2 (\vec{\upsilon}\_i \times \vec{\upsilon}\_k) = \vec{\upsilon}\_i \times \vec{\upsilon}\_k.$$

Since the direction of the static current obviously has no effect on the cross-product, the result for the yaw angle also remains unchanged. Flipping the signs of the field components in Equation (10) for the pitch angle leads to:

$$\frac{-\upsilon\_{kx} - \upsilon\_{ky}A}{-\upsilon\_{kz}\sqrt{A^2 + 1}} = \frac{\upsilon\_{kx} + \upsilon\_{ky}A}{\upsilon\_{kz}\sqrt{A^2 + 1}}, \quad A = \tan\alpha.$$

The direction of the static current obviously also has no influence on the pitch angle.

Equation (14) for the roll angle *γ* uses the measurement from only one sensor as the input. Thus, it is also independent of the sign of *I*. In summary, the sign of *I* does not affect the orientation as long as the measurements are referenced to the same time interval of the corresponding signal.

#### **4. Simulation Results**

In this section, the presented equations for determining the UAV position and orientation are tested in the context of a software-in-the-loop (SIL) simulation.

#### *4.1. SIL Design and Test Procedure*

The SIL simulation was implemented in MATLAB and cyclically ran through a specified flight trajectory created with the Gazebo tool. Each waypoint of the trajectory is represented by a vector (*yD*, *zD*, *α*ˆ, *β*ˆ, *γ*ˆ)*T*, which describes the position and orientation of the local UAV frame *Lu* with respect to the global frame {*G*}. From this, the SIL computes the positions of the virtual sensors *i* and the corresponding magnetic measurements *vi* with respect to the local sensor frame *Li* in each run. For this purpose, Equations (4) and (5) were used, where the rotation matrix was used according to the convention in Gazebo, i.e., *R*ˆ *LU G* = *RT x* (*γ*ˆ)*RT y* (*β*ˆ)*RT z* (*α*ˆ). This corresponds to the transposed version of the definition used in this paper, i.e., *R*ˆ *LU G* = (*RLu G* )*T*. In order for the angles *α*, *β*, and *γ* obtained in the test run to be compared later with the Gazebo angles *α*ˆ, *β*ˆ, and *γ*ˆ, they must first be converted at the end of each run. This is done automatically by the SIL. To calculate the measured values *vi*, the SIL assumes that there is a transmission line system with the arrangement described in Section 2.1; i.e., the cables are in the *xy*-plane at −*y*0 and *y*0. The SIL then passes the computed *vi* as the input values to the test procedure, which verifies the entire method described in this publication. The test procedure uses the constraints defined in Section 3.3 for navigation as a priori knowledge. It is called cyclically from the SIL until all waypoints on the flight trajectory have been passed and is divided into the following steps:

$$
\delta L\_k = \sum\_{i=1}^{N} (\vec{p}\_{i,\varepsilon} - \vec{p}\_{i,k})^T (\vec{p}\_{i,\varepsilon} - \vec{p}\_{i,k}),
\tag{32}
$$

where *N* is the number of nodes, respectively sensors. Choose the polygons of the two smallest sums. Finally, choose one out of these two whose perimeter most closely matches that of the original polygon.

10. Calculate the UAV position from the found polygon for each node using Equation (4), i.e.,

$$
\vec{p}\_{D,i} = \vec{p}\_{i,k} - \mathcal{R}\_G^{L\_u}(\vec{u}, \vec{\beta}, \vec{\gamma}) \,\vec{p}\_{i0}.
$$

Take the average of all *pD*,*i* as the UAV position *pD*.

If the entire trajectory has been run through by the SIL, the specified and actual position are output graphically. The same goes for the angles. The software offers several configuration parameters. These are:

#### *4.2. Test Parameters*

The test cases and the associated results are presented below. If not otherwise stated, the following configuration is used:

$$
\vec{p}\_1 = \begin{pmatrix} 0.05m\_\prime 0.2m\_\prime - 0.05 \end{pmatrix}^T, \quad \vec{p}\_2 = \begin{pmatrix} 0.05m\_\prime - 0.2m\_\prime - 0.05 \end{pmatrix}^T.
$$

• Two sensors rotated by 90◦. Their positions in *Lu* are:

$$
\vec{p}\_3 = (0.05m, 0.2m, -0.25)^T, \quad \vec{p}\_4 = (0.05m, -0.2m, -0.25)^T.
$$

In the previous publication [25], it was shown that the electric motor of a small 1:5 scale electric toy car can produce stray magnetic fields of several μT; especially when the vehicle accelerates. The reason for this is that inadequate magnetic shielding measures are often taken in technical installations. The supply lines of electronic components with high power requirements can therefore emit strong magnetic interference fields. Since the UAV was also equipped with several electric motors, the test cases were performed with different standard deviations *σ* of the white noise. A value of 1 μ*T* was assumed for the case of good shielding and 8 μ*T* for poorly shielded electric motors. These were realistic values based on the measurements in [25] and took into account that the UAV may have multiple electric motors.

#### *4.3. Test Results*

Figure 10 shows the localization results using the derived equations and the proposed polygon detection method.

**Figure 10.** Localization quality, good shielding: *σ* = 1 μT.

The flight started at the arrow in the indicated direction, with the UAV flying one lap under the transmission lines. As was to be expected, good shielding allowed more precise navigation. The quality of localization was obviously very high with a magnetic noise standard deviation of *σ* = 1 μT. The mean deviation from the specified flight path was 21 mm, with a maximum error of 54 mm. The standard deviation of the roll angle error was 0.4◦, and its mean was very close to zero, around 0.025◦. Only the roll angle is shown here, as it had the largest values since the UAV moved almost exclusively in the *yz*-plane in the simulation.

Figure 11 shows the localization results for the case where poorly shielded electrical drives in the UAV produce magnetic noise with a standard deviation of *σ* = 8 μT.

**Figure 11.** Localization quality, poor shielding: *σ* = 8 μT.

Despite the high magnetic interference, the navigation worked very well near the transmission lines due to the high signal-to-noise ratio (SNR). The deviations increased with greater distance from the transmission lines. Figure 12 shows the deviation from the specified flight path and the specified roll angle for each waypoint.

**Figure 12.** Path and roll angle errors, poor shielding: *σ* = 8 μT.

The mean deviation from the specified flight path was 41 mm, where the maximum error was 221 mm. The mean deviation from the specified roll angle was −0.22◦, where the maximum error was 12.3◦. The standard deviation of the roll angle error was 2.96◦.

Figure 13 shows the localization results for the same noise standard deviation of *σ* = 8 μT when only three sensors were used, i.e., Sensor 4 at *p*4 was excluded.

**Figure 13.** Localization quality, three sensors, *σ* = 8 μT.

The mean deviation from the specified flight path was 91 mm. In areas far from the transmission lines, as in this case, very strong outliers could occur occasionally, so that the maximum error was now about 2260 mm. The outliers occurred very frequently compared to the case with four sensors, as demonstrated by repeated experiments. These were not due to an incorrect calculation of the UAV position, but to an incorrect detection of the polygons. These were extremely distorted by the strong magnetic interference and could therefore hardly be distinguished from each other. This was to be expected, because the more sensors were used, the higher the noise immunity was. In practice, the UAV would carry out aerial manipulations in the vicinity of the transmission lines and would therefore be in an area with a high SNR. Therefore, it may also be possible to use only three sensors, as shown in Figure 14. The UAV started its flight near the transmission lines and then flew between the lines. Only when it then moved away did strong outliers occur. Repeated tests always led to the same result.

**Figure 14.** Localization quality close to the transmission lines and far away, three sensors, *σ* = 8 μT.

In the tests, a current of 700 A was assumed. It should be mentioned that the quality of the localization depended on the current strength in the transmission lines. This in turn depended on the per capita energy consumption of the residents, which may vary from city to city. Thus, for completeness, Figure 15 shows additionally the same case as Figure 10 with a low transmission line current of 100 A.

**Figure 15.** Localization quality in the case of low transmission line current of 100 A, four sensors, *σ* = 1 μT.

#### **5. Discussion**

As mentioned in Section 3.2, in a real transmission line system, it is difficult to detect whether the UAV is located below or above the cables without additional detection mechanisms. This is only possible in Figure 14 because the SIL calculated with a static equivalent current and was only used in the context of the tests.

In this paper, transmission lines were modeled as perfectly straight cables, which may not be the case in real-world circumstances. Transmission lines, for example, hang downward, thus forming more of a curved line. Even though the radius of this curve is very large, the influence on navigation accuracy should be investigated in a further work. In UAV indoor navigation applications, this problem does not occur because the cables are shorter and can be routed to fit the straight line model.

Furthermore, the magnetic field contributions of the other two phases of the threephase system were neglected, since they are typically at least 2 m away from the phase under consideration. This should also be investigated in more detail in a further study in the context of transmission line navigation.

Finally, the magnetometers should be as far away as possible from ferromagnetic materials that may be in the UAV; otherwise, they may affect the sensor readings.

#### **6. Conclusions**

In this work, it was demonstrated that UAV navigation is possible in the magnetic field of two long parallel conductors. In this context, equations were derived to be able to calculate the UAV orientation analytically. It was shown that only three magnetometers are needed for this purpose and that neither the current, nor the magnetic field equations

need to be known. Furthermore, it was shown that combining the magnetic field equation with the magnetic energy equation leads to analytical solutions for the calculation of the sensor position. This results in 16 possible positions for a single sensor, from which the correct one must be determined. For this purpose, theorems were derived, which can be used to exclude 14 of the 16. The last two remaining positions are not distinguishable. However, it was shown that considering multiple sensors allows a whole set of correct positions to be selected. Here, the positions were combined into a number of polygons in the *yz*-plane, from which the one most similar to the original polygon was selected. This can be determined from the known orientation and the user-defined sensor positions in the local UAV frame. This ultimately led to a unique determination of the correct sensor positions and hence the UAV position. All equations, as well as the polygon selection method were successfully verified for different noise.

**Author Contributions:** Investigation, D.M.; methodology, D.M.; software, D.M.; validation, D.M.; writing—original draft preparation, D.M.; writing—review and editing, Z.K., S.B.; visualization, D.M.; supervision, S.B., Z.K.; project administration, S.B. All authors read and agreed to the published version of the manuscript.

**Funding:** The project AERIAL COgnitive integrated multi-task Robotic system with Extended operation range and safety (Aerial-Core) received funding by the European Union's Horizon 2020 research and innovation program under grant agreement No. 871479.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** I would like to thank my colleague Filip Zori´c who generated the test flight data in Gazebo and thus supported the tests.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**