Virginie Griseri

Fulbert Baudoin

Jizzini, Hadi, Yussef, jazoo, Farhat, Jihan, Jinan, Maha, Hsein, Zeinab, Sakka, Kolkas, Kassem, AbuZain, Daher Jana Hassan

Florian Bugarin and Stephane Segonds for their endless help, advices

Introduction

Electrical insulators made of polymeric materials are widely used in various applications, from electronic components to electrical machines. Among the electrical insulators with a complex formulation, this work will focus on Low-Density Polyethylene (LDPE). This work addresses the transport of electrical energy in high-voltage cables under DC stress. This electrical stress induces internal charges called space charge, whose accumulation is directly linked to the degradation of the dielectric, leading to the electrical insulator breakdown. The major difficulty arises from the fact that phenomena describing the dielectric aging is not fully-understood. Therefore, Additional studies are required to understand this charge accumulation and the mechanisms by which the insulator is affected to develop strategies for material improvement and the formulation of behavioral laws under stress.

There exist two main approaches that could be useful for these studies; they are known as experimental and numerical modeling approaches. The originality of this study lies in the intention to enhance the numerical modeling approach by integrating it with the experimental one using optimization tools, which can provide many beneficial information.

Most of the actual physical concepts used right now to describe charge transport in solid dielectrics have been studied for more than 30 years. Two physical models have been developed in our laboratory to describe the Bipolar Charge Transport (BCT) in LDPE under DC stress. These models are known as fluid models, and they are defined as follows:

• Model 1: This model considers two levels of traps (shallow and deep traps), as well as the most essential phenomena taking place inside dielectrics, such as injection, hopping mobility, and the recombination process of positive and negative carriers. Two kinds of carriers are considered in this model, mobile or trapped. These carriers are generated by a modified-Schottky injection at the metal-dielectric interface and are extracted without a barrier. The model is based on Poisson's equation and the conservation law of charges.

• Model 2: This model is similar to Model 1; however, the main feature of this model is that trap depths are represented by an exponential distribution, instead of considering only two levels of traps. Besides, all charges are supposed to be trapped with no mobile charges.

Since the scientific approach addressed in this study is innovative in our research community, it is preferred to start with a model considering the physical phenomena that are irrelevant to each other and capable of describing the charge transport in LDPE under DC stress. For this reason, we adopted Model 1 for this study instead of using Model 2. This model requires a large number of unknown parameters that represent the physical phenomena taking place in dielectrics, such as injection barrier, mobility, trapping, and de-trapping coefficients. Most of these parameters cannot be predicted, observed, or estimated by independent experiments. For this reason, an optimization algorithm is used to optimize the BCT model to provide a realistic representation of the physical, electrical, and chemical phenomena of dielectrics under stress. Our study employed this algorithm to find the best model parameters by minimizing the sum of the squares of the deviations between the experimental and simulated data. Experimental data that have been used are the net density of charges measured by the Pulsed Electro-Acoustic method (PEA) along with what is known as charging current measurements. On the other side, the simulated density of charge and current could be obtained by the BCT model developed for LDPE under DC stress. All experiments and samples preparations were performed at LAPLACE Laboratory. This thesis is composed of five chapters, as follows:

Chapter 1 presents the state of art of this study, which highlights the essential physical phenomena used to describe the charge transport in LDPE. The first part of this chapter presents the chemical and physical properties of polymers, particularly LDPE. Besides, the charge transport and charge generation phenomena in LDPE will be addressed in detail. Furthermore, the mathematical model used in this study is proposed to describe the charge transport and electrical behavior in LDPE, which addresses most of the microscopic and macroscopic phenomena during the polarization and depolarization processes. The second part of this chapter shows the contributions made by several researchers over the previous 30 years to the BCT model. In conclusion, the main objective of the thesis will be presented by specifying the issues we are trying to solve. Numerical optimization is not widely used in our community of research. For this reason, Chapter 2 reviews the basic concepts of optimization that should be well understood in order to apply optimization to a specific problem, such as, global and local optimality conditions. The general form of an optimization problem can be classified according to the nature of the objective function and constraints. Thus, many optimization models depending on the type of the cost functions and constraints will be presented and discussed in order to know which one represents the problem we are attempting to solve.

The main goal of Chapter 3 is to reveal the most suitable optimization algorithm for solving our problem. The first part provides a detailed investigation of optimality requirements that demonstrates the strong and weak points of the most significant algorithms that could be applied for optimizing the BCT model. Five optimization algorithms will be studied for estimating the unknown set of parameters related to the BCT model. All the algorithms will be compared in detail by examining their robustness, efficiency, and accuracy when applied to the BCT model. Then the best algorithm will be chosen to be used in the following sections.

Chapter 4 proposes an original method based on a dual approach Simulation/Experiment that helps to optimize and critically analyze the BCT model (i.e., highlighting the weaknesses and strengths). This approach aims at constructing a multi-objective cost function that integrates both current and charge experiments with varying electrical fields in a single cost function. Then, this cost function will be minimized using the algorithm chosen in Chapter 3. Besides, the sample conditioning and the experimental setups utilized for this study will be discussed. Since this study aims to combine both current and charge experiments in a single cost function, the samples were carefully manufactured with specific characteristics to suit both experiments. Finally, the results will be discussed and analyzed in depth.

Finally, Chapter 5 focuses on studying and analyzing the charge injection phenomena taking place at the electrodes using optimization tools. Several types of surface electrodes will be considered in this chapter, that are: gold, aluminum, copper, and semi-conductor. This chapter intends to compare the nature of different electrodes by classifying the amount of injected charges through each one. However, it is quite difficult to obtain an accurate value of the injected charges during voltage application using the PEA method. Yet, it is possible to roughly estimate and compare the quantity of injected charges at short time, and then identify the influence of the surface electrode type. As in Chapter 4, optimization tools will be used to find the unknown parameters representing the injection barrier height that can fit experiments with model simulations. Finally, the experimental measurements will be compared with simulated results obtained by the BCT model to highlight the reliability of the presented method.

Chapter I

I.State of Art

I.1 Introduction

The vast majority of the actual physical concepts used to depict charge transport and breakdown in solid dielectrics, have been studied for more than 30 years. These studies provided a wealth of information regarding amorphous semiconductor [1], charge localization and transport [2], charge states in polymers [3], hopping conduction, space charge limited current, interface states, etc. Some great review articles depict the basis of our present information in this area [4]- [10]. T. J.

Lewis [4], [8] discussed charge transport, charge injection, and breakdown in polymeric insulators, especially in polyethylene. H. J. Wintle [5], [7] studied the basic physics of insulators, such as, charge motion and trapping phenomena. T. Mizutani [9] and S. Boggs [10] addressed the high field phenomena in insulating polymers.

Despite this physical mechanism background, studies on charge transport modeling have been limited over the past 30 years. Most of the time, they have considered analytical models in simple case studies or under prohibitive conditions rather than real conditions of materials in electrical systems. The trouble for advancing such models has been the absence of essential information, such as carriers traps, charge mobility and their local density, etc. Besides, information on the microscopic process controlling charge transport was inaccessible, since all the experimental measurements depended on strategies said to be "integral", which acquired the information in time and/or in space (for instance: potential surface measurement of a dielectric resulting from an internal distribution of charges). Two advancements have contributed to change the situation:

• the first is due to the enormous increase of the computers calculating power and to the improvement of computational techniques [5];

• while the second concerns the advancement of the experimental techniques which today permit the measurement of the internal distribution of the space charge as a function of time.

The birth and the development of these measurement techniques indicate a significant advance in the dielectric area. Several methods have arisen to measure the space charge on polymers and have provided new information, especially on Polyethylene (PE) [11]- [16]. These methods provide some new information, such as the carriers' polarity, mobility of charge carriers, or the trap's depth Chapter I 20 using different electric fields and temperatures. Most information on recombination processes can be obtained by combining space charge calculation techniques with various spectroscopic techniques, showing whether the transport is bipolar or unipolar [17], [18]. Furthermore, other factors that have impressively changed the situation and provided a boost to the improvement of models describing the charge transport phenomena within insulating polymers under stress are:

• The evolution towards a compact system in power electronics and, more generally, in all the domains of electrical engineering, leads to an increase in power density, with two consequences for polymeric insulations: first, these materials are consistently exposed to strong stresses, often close to their limits; second, manufacturers are developing new instruments capable of making best estimations of the intrinsic limits of the insulator used to ensure the security of their usage;

• The demand for safer electrical systems is principally due to their use in sensitive applications (large-scale electrical networks, embedded networks, complex systems, etc.).

That is why new predictive models that anticipate the effects of chemical and/or physical defects should be formulated;

• The development of new materials for electrical insulation, such as, composites with inorganic charges on a micro-or nano-metric scale [19], with chemical functionalization (by addition or grafting) of the matrix or the charges. The models must be appropriately formulated to fit these results.

The aim of this chapter is not to evaluate previously published modeling works or to introduce new dielectric physical concepts. We intend to show that by considering the accumulated knowledge in transport physics, there is a major potential to improve simulation practices to model the behavior of polymer insulators under electrical stress. Hence, this work tries to provide a realistic representation of physical, electrical, and chemical phenomena at the microscopic scale that are not properly handled in macroscopic models. The term "modeling" is used here in its general context; "some techniques will be based on microscopic principles, while others may be based on macroscopic models, depending on the nature of the problem''. These models provide a wide range of data that can be used for a variety of purposes, such as further understanding the mechanics of insulating materials or developing new scientific and industrial applications.

Chapter I 21

In the first section of the chapter, we will discuss the chemical and physical properties of polymers, particularly Low-Density Polyethylene (LDPE). The basics of molecular properties, such as defects, trap levels, and carrier's generation will be described. In the second part, we will describe the space charge in polymers. This part will address the charge transport and charge generation phenomena taking place in LDPE. In the third section, a mathematical model is proposed to describe the charge transport and electrical behavior in LDPE. This model addresses most of the microscopic and macroscopic phenomena taking place in LDPE during the polarization and depolarization processes, such as injection, recombination, trapping, and detrapping phenomena. Also, the evolution process of this model is discussed by showing the contributions made by several researchers over the previous 30 years. Finally, the main objective of the PhD is presented by specifying the issues we are trying to solve and the solutions that will be addressed in the following sections.

Among insulating materials with a large bandgap, polyethylene (PE) stands out for its simple chemical structure (repetitive unit CH 2 ), its chemical inertness, and its numerous electrical engineering applications. This material has been most often considered for several fundamental studies (molecular descriptions, macroscopic modeling), making it by far the most studied material.

I.2 Chemical and physical properties of polymers I.2.1 Polymers

A polymer is a long organic chain consisting of smaller repetitive simple molecules, called monomeric molecules. The electrical properties of these macromolecules depend on the chemical nature of the chain units, their design, and eventually their operating conditions. The question is to determine the relationship between the material's electrical properties and the microscopic structure. At the macroscopic scale, a polymer is considered a homogenous substance. Still, at a smaller scale, its structure is significantly more complicated due to the number of constituent atoms, the arrangement of these atoms between them, and implementation parameters. The chemical and physical structure of a polymer will then be defined. Then, in terms of structural features and transport processes, we'll concentrate more particularly on the properties of polyethylene, our research material.

Chapter I

22

I.2.1.1 Chemical structure

A polymer is generated by linking the repeated units together through covalent connections (which come from forming a bond between the monomers). The polymer backbone is made up of carbon chains. The monomers may be inserted in a straight line or a branching pattern during the polymerization process. These two types of configurations are called linear and branched. Such structures are associated with the bulk of the initial monomer and/or chemical products, making it possible to connect chains laterally to the macromolecule skeleton. The synthesis of a polymer is achieved in various ways: anionic polymerization, cationic, radical, emulsion polymerization. The common point among these syntheses is the presence of initiators and/or catalysts, which are added to increase polymerization speed. Unfortunately, the "pure" polymers are fragile. They can suffer degradation, which leads to a change in their chemical structure and physical properties under different sources (heat, moisture, light, etc.). Therefore, they cannot be used as they are. During their development, the manufacturers add various substances which modify their properties:

• Anti-oxidants, which prevent polymers from reacting in the presence of oxygen in the air • Cross-linking mechanisms, which after a chemical reaction can bind to the polymer chains to form inter-molecular bridges. These agents allow the formation of a macromolecular network, which gives the material a better mechanical and thermal stability.

• Other additives can also be added, such as plasticizers, dyes etc., depending on the properties required.

In the process of polymerization, all those additional additives can react with other molecules, and sample substrates can contain their degradation products.

It is difficult to classify or measure all impurities in the substructure, but there is a particular discrepancy in these products. We call this disorder of chemical type because residues from the polymer's synthesis primarily cause it. We will later see the role these impurities can play in the charge transport in these materials.

I.2.1.2 Physical structure

By physical structure, we mean the architecture of the chain during the elaboration of the polymer.

The simplest arrangement is the arrangement of the atoms or groups of atoms with each other.

These can be arranged by rotating around the covalent bonds to adopt the conformation that requires the least energy.

Chapter I

23

The orderly character of molecules or atoms over large distances is referred to as crystallinity.

There are often defects, such as ramifications, and the polymer only partially crystallizes. It is referred to as semi-crystalline. It has two phases: a crystalline, well-ordered phase and an amorphous phase. If no structure exists, the polymer is amorphous.

As discussed in the previous paragraph, chemical impurities deform the crystallographic mesh, modifying the physical structure of the polymer. These physical structure distortions result in an imbalance within the polymer structure, which we can refer to as a physical or conformational disorder.

I.2.2 Defects in polymers

Low density polyethylene is semi-crystalline, consisting of a crystalline part organized in the form of lamellae of approximately 10 nm in size, and an amorphous part. The lamellar molecular arrangement is shown in Figure 1a. The amorphous part is due to the conformational and chemical disorder, introduced among others by the unsaturated species as shown in Figure 1b. b) examples of unsaturated species at the origin of the chemical disorder. [24] LDPE materials contain different structural irregularities and defects, such as branch points, conformational defects, or cross-links. Such defects may arise due to the faulty chemical process of linking monomer units, yielding an undesired linkage and breakage of chemical regularity [25].

Another chemical defect is branching, which occurs when a linear polymer chain splits into two separate chains that occur at different rates. Furthermore, cross-links produce unique mechanical properties, e.g., rubber elasticity is due to cross-linking [START_REF] Zerbi | Transport of matter upon annealing in crystalline polymethylene systems: a calorimetric and spectroscopic study[END_REF]. Defects can also occur in the conformational coiling of polymer chains [25]. Conformational defects arise due to kinetic and energetic influences in the crystallization of an otherwise chemical and stereo-chemically regular chain. Conformational errors may also be strongly caused by chemical or stereo-chemical defects;

induced conformational defects are also named.

Chapter I

24

Defects in polymers have been studied using several experimental techniques, including Scanning Electron Microscopy (SEM), image analyzer, and Transmission Electron Microscopy (TEM) [START_REF] Malhotra | Recent studies of heterocyclic and aromatic conducting polymers[END_REF].

With SEM, it has been possible to determine the voids, grain-size, and free-volume in polymeric materials. X-ray diffraction measurements were helpful in imaging vacancies, clusters, and dislocations in the crystalline regions of polymeric materials.

The temperature dependence of density, electrical conductivity, and the dielectric constant, as well as the frequency dependence of each of these parameters, can give information on the dynamic nature of the defects. Additionally, there have been problems associated with isolating situations in which specific structural defects have a clearly defined role in controlling macroscopic properties. In the same manner, the study of defects in polymers is a reasonable source that provides a relation between microscopic and macroscopic structural properties that are the fundamental theme of all material science.

Obviously, it is quite difficult to obtain a clear observation of the chemical and physical defects in polymers. For this reason, in this work, one of the essential targets is to better understand the defects by using a dual approach (experimental/modeling).

I.2.3 Trap sites

One of the major goals of this research is to investigate the defects in conducting polymers. The chemical and physical structures of the material might be used to analyze these defects. Trapping levels are commonly used to characterize defects inside a material. Only two trapping levels will be addressed in our model: shallow traps, which represent the material's physical defects, and deep traps, which represent the material's chemical defects.

I.2.3.1 Shallow traps

In general, the shallow trap states represent the physical defects and structural inconsistencies inside the material [START_REF] Teyssedre | Charge transport modeling in insulating polymers: from molecular to macroscopic scale[END_REF], [START_REF] Meunier | Molecular modeling of electron traps in polymer insulators: Chemical defects and impurities[END_REF], which produce a reduction in the carrier mobility. Shallow traps feature energetic depths in the range of 0.05-0.75 eV (Figure 2) [START_REF] Meunier | Molecular modeling of electron trapping in polymer insulators[END_REF], [START_REF] Cubero | Electronic states for excess electrons in polyethylene compared to long-chain alkanes[END_REF], (measured from the original band edge) and have high densities in the range of 1026-1027 m -3 [START_REF] Teyssedre | Charge transport modeling in insulating polymers: from molecular to macroscopic scale[END_REF], [START_REF] Meunier | Molecular modeling of electron traps in polymer insulators: Chemical defects and impurities[END_REF], [START_REF] Sato | Quantum chemical calculation of hole transport properties in crystalline polyethylene[END_REF]. The carriers residence time in the shallow traps is on the order of 10 -11 to 10 -13 s, which tends to be consistent with the fact that these sites support transport but do not allow stability of the charges for the long-term that are detected using traditional techniques. For instance, based on the authors, we find several values for the average inter-trap value [START_REF] Sato | Quantum chemical calculation of hole transport properties in crystalline polyethylene[END_REF]- [START_REF] Roy | Modelling space charge in a cable geometry[END_REF].

Chapter I

25

In the material, the major mode of spatial charge transport takes place in the shallow trap states.

Yet, due to their low energetic barrier and minor spacing, transport in between them occurs at a particularly high speed, so that carriers within shallow traps cannot be measured with PEA nor productively represented in the simulation model.

I.2.3.2 Deep traps

Deep states are physio-chemical defects inside the material that allow the material to display slow charge decay [START_REF] Meunier | Molecular modeling of electron traps in polymer insulators: Chemical defects and impurities[END_REF], [START_REF] Meunier | Molecular modeling of electron trapping in polymer insulators[END_REF]. Deep traps are known to have energetic depths, 𝑤 𝑡𝑟 in the range of 0.9-2 eV (Figure 2) with a density in the range of 10 20 -10 21 𝑚 -3 [START_REF] Teyssedre | Charge transport modeling in insulating polymers: from molecular to macroscopic scale[END_REF], [START_REF] Boufayed | Models of bipolar charge transport in polyethylene[END_REF]. Because of their low density, inter-deep trap conduction without the interaction of inter-level and shallow traps is highly improbable. In our manuscript, deeply trapped carriers are referred to as "trapped".

The residence time of the charges in the deep traps is much higher than those encountered for shallow traps, in the order of 4.6 × 10 -3 s to 5.1 × 10 5 s [START_REF] Gao | Charge transport in low density polyethylene based micro/nanocomposite with improved thermal conductivity[END_REF], and charges can be stabilized by deep traps for virtually infinite durations.

I.2.4 Polyethylene (LDPE)

The polymer we have chosen to study is polyethylene. It is a member of the vinyl polymer family and is derived from the polymerization of ethylene, or ethene (Figure 3). The repetition of the basic unit makes it possible to obtain polyethylene of chemical formula [-C2H4-]n (Figure 4), where n represents the degree of polymerization of the molecule.

Chapter I Its electrical and mechanical qualities (Table 1) have made it a preferred material for high voltage cable insulation. Polyethylene (PE) has good insulating properties such as resistivity, dielectric strength, and dielectric losses. Hence, it has been widely used in a large variety of applications in the electrical and electronic industry. As the primary insulation of a power cable, the PE material could be subjected to a large amount of heat when high power is transmitted through the cable line. Furthermore, due to the low thermal conductivity of the PE matrix, the phenomenon of heat accumulation is quite severe, which could accelerate the thermal degradation of cable insulation and lead to dielectric breakdown [20]- [23], posing a threat to the power system's safe operation. However, LDPE is simpler to process, and its cost is much lower. Therefore, the LDPE was chosen to be studied in our research.

I.3 Space charge generation and transport

Space charge is a concept in which excess electric charge is treated as a continuum of charge distributed over a region of space (either a volume or an area) rather than distinct point-like charges. This model typically applies when charge carriers have been emitted from some region of a solid and if they are sufficiently spread, or the charged atoms or molecules left behind in the solid can form a space charge region. The space charge phenomena could be described by addressing two major topics: how carriers are generated and how they are transported. The following sections describe the two processes in detail.

I.3.1 Charge generation

The generation of charge mechanisms are important to the successful applications of dielectric materials because it is challenging to access phenomena occurring at a microscopic scale. There are three kinds of generation mechanisms in dielectric films, that are:

• irradiation

• ionization

• injection conduction mechanism

Electron-beam irradiation is an alternative way to generate charges in insulating materials at controlled positions and quantities in order to monitor their behavior regarding transport phenomena under the space charge induced electric field or external field applied. Similarly, charges could be generated when an electric charge is ionized. Ionization is the process by which an atom or a molecule acquires a negative or positive charge by gaining or losing electrons, often in conjunction with other chemical changes. The mechanisms of charge generation at the interfaces are known as "injection conduction." Numerous studies conducted on this subject have provided a partial understanding of these interface-controlled processes, particularly for metal-insulator type interfaces. In our case, the charges are assumed to be produced only by injection at both electrodes, with no consideration given to irradiation or ionization processes.

The injection conduction mechanisms depend on the electrical properties at the electrode-dielectric contact. The barrier height at the electrode dielectric interface is the most important parameter in this type of conduction mechanism. The injection conduction mechanisms include:

• Schottky or thermionic emission

• Fowler-Nordheim tunneling,

• thermionic-field emission.

The current caused by thermionic emission is strongly dependent on the temperature, whereas the tunneling current is almost temperature independent. Aside from the barrier height at the electrodedielectric interface, the effective mass of the conduction carriers in dielectrics is also a critical feature for the injection conduction mechanisms.

I.3.1.1 Schottky or thermionic emission.

Schottky emission is a conduction mechanism at the electrodes [START_REF] Dissado | Elemental strain and trapped space charge in thermoelectrical aging of insulating materials. Part 1: Elemental strain under thermo-electrical-mechanical stress[END_REF], [START_REF] O'dwyer | The theory of electrical conduction and breakdown in solid dielectrics[END_REF], where the electrons in the metal will overcome the energy barrier at the metal-dielectric contact, if they can obtain enough energy from thermal activation to pass into the dielectric. Thermionic emission is one of the most frequently observed conduction mechanisms in dielectric films, especially at relatively high temperatures. Figure 5 shows the metal-insulator energy band diagram when the metal electrode is under negative bias with respect to the dielectric and the semiconductor substrate. The image force can reduce the energy barrier height at the metal-dielectric contact. The barrier-lowering effect caused by the image force is defined as the Schottky effect. Such conduction mechanism due to electron emission from the metal to the dielectric is called thermionic emission or Schottky emission. The expression of Schottky emission is

𝐽 𝑖𝑛𝑗 = 𝐴 • 𝑇 2 • 𝑒𝑥𝑝 (- 𝑞 • 𝑤 𝐵 𝑘 𝐵 • 𝑇 ) [𝑒𝑥𝑝 ( 𝑞 𝑘 𝐵 • 𝑇 √ 𝑞 ⋅ 𝐸 4𝜋𝜀 0 𝜀 𝑟 )] (1)
where, 𝐽 𝑖𝑛𝑗 is the current density; D is the thickness of the material; T is the absolute temperature, 𝑞 = 1.6 × 10 -19 𝐶 is the elementary charge, 𝑤 𝐵 is the injection barrier, respectively, 𝑘 𝐵 stands for Boltzmann constant, 𝐸 is the electric field, 𝜀 0 = 8.85 × 10 -12 𝐹. 𝑚 -1 is the permittivity of vacuum and 𝜀 𝑟 = 2.3 is the material relative permittivity.

Chapter I

𝐴 = 4𝜋𝑞𝑘 2 𝑚 * ℎ 3 = 120𝑚 * 𝑚 0 = 1.2 × 10 6 𝐴𝑚 -2 𝐾 -2
Where 𝐴 is the effective Richardson constant, 𝑚 0 is the free electron mass, 𝑚 * is the effective electron mass in dielectric, ℎ is the Planck's constant.

I.3.1.2 Fowler-Nordheim tunneling.

According to classical physics, the electrons will be reflected if their energy is less than that of the potential barrier. However, quantum mechanics predicts that when the potential barrier is thin enough (<10 nm), the electron wave function will pass through. Hence, due to the tunneling effect, the probability of electrons being on the opposite side of the potential barrier is not zero. Figure 6 represents the schematic energy band diagram of Fowler-Nordheim (F-N) tunneling. F-N tunneling arises when the applied electric field is strong enough to allow the electron wave function to pass through the potential barrier and into the dielectric's conduction band. The F-N tunneling current is expressed as:

𝐽 𝑖𝑛𝑗 = 𝑞 3 𝐸 2 8𝜋ℎ𝑞𝑤 𝐵 exp [- 8𝜋(2𝑞𝑚 𝑇 * ) 1 2 3ℎ𝐸 𝑤 𝐵 3/2 ] (2)
where 𝑚 𝑇 * is the tunneling effective mass in dielectric; the other notations are the same as defined before. The tunneling current may be computed by measuring the current-voltage (I-V) properties of the devices at very low temperatures. At such a low temperature, the thermionic emission becomes negligible and the tunneling injection is dominant. emission is illustrated in Figure 7. The current density due to thermionic-field emission can be expressed as [START_REF] Hamann | Electrical conduction mechanisms in solids[END_REF]:

𝐽 = 𝑞 2 √𝑚 (𝑘𝑇) 1 2 𝐸 8ℎ 2 𝜋 5/2 exp(- 𝑞𝑤 𝐵 𝑘𝑇 ) exp [ ℎ 2 𝑞 2 𝐸 2 24𝑚(𝑘𝑇) 3 ]
(3) Poole-Frenkel (P-F) emission entails a mechanism that is very similar to Schottky emission; in particular, the thermal excitation of electrons may emit from trap states into the conduction band of the dielectric. Also, the P-F emission is sometimes entitled as the internal Schottky emission.

Considering an electron in a trapped state, an applied electric field across the dielectric film can Chapter I 32 reduce the electron's Coulomb potential energy. The reduction in potential energy may enhance the probability of an electron being thermally excited out of the trap and into the dielectric's conduction band. The schematic energy band diagram of P-F emission is shown in Figure 9. For a Coulombic attraction potential among electrons and traps, the current density due to the P-F emission is

𝐽 = 𝑞𝜇𝑁 𝐶 𝐸 exp [ -𝑞 (𝜙 𝑇 -√ 𝑞𝐸 𝜋𝜀 𝑟 𝜀 0 ) 𝑘 𝑏 𝑇 ] ( 4
)
Where 𝜇 is the electronic drift mobility, 𝑁 𝐶 is the density of states in the conduction band, 𝜙 𝑇 is the trap energy level, and the other notations are the same as defined before. Because P-F emission is caused by thermal activation in the presence of an electric field, this conduction mechanism is frequently observed at high temperatures and strong electric fields.

I.3.2.2 Hopping Conduction.

Hopping conduction is made by thermal excitation above a potential barrier, by quantum tunneling effect of trapped electrons "hopping" from one trap site to another in dielectric films. Figure 10 shows the schematic energy band diagram of hopping conduction. The expression of hopping conduction is [START_REF]Handbook of nanoceramics and their based nanodevices[END_REF]- [START_REF] Mott | Electronic processes in non-crystalline materials[END_REF]:

𝐽 = 𝑞𝑎𝑛𝑣 exp ( 𝑞𝑎𝐸 𝑘 𝑏 𝑇 ) exp (- 𝑞𝜙 𝑇 𝑘 𝑏 𝑇 ) ( 5
)
where 𝑎 is the mean hopping distance (i.e., the mean spacing between trap sites), 𝑣 is the frequency of thermal vibration of electrons at trap sites, and 𝜙 𝑇 is the energy level from the trap states to the bottom of conduction band (𝐸 𝐶 ); 𝑛 is the electron density in the conduction band of the dielectric, the other terms are as defined above. The P-F emission corresponds to the thermionic effect and the hopping conduction corresponds to the tunnel effect. The thermionic mechanism allows carriers to cross the trap barrier in P-F emission. However, in hopping conduction, the carrier energy is lower than the maximum energy of the potential barrier between two trapping sites. In this situation, the carriers can still use the tunnel mechanism to move.

I.3.2.3 ohmic conduction

Ohmic conduction is produced by the movement of mobile electrons and holes in the conduction band and valence band, respectively. In ohmic conduction, a linear relationship occurs between the electric field and current density. Figure 11 shows a schematic energy band diagram of the Ohmic conduction. Despite the fact that dielectrics have a high energy band gap by definition, thermal excitation will generate a small number of carriers. For example, electrons from the valence band or the impurity level may be excited to the conduction band. The carrier numbers will be extremely low but they are not zero. The current density of ohmic conduction can be defined as:

𝑗 = 𝜎𝐸 = 𝑞𝜇𝐸𝑛 (6)
Chapter I

34

The magnitude of this current is very small. This current mechanism may be observed when there is no significant contribution from other conduction mechanisms of current transport in dielectrics [START_REF]Handbook of nanoceramics and their based nanodevices[END_REF].

I.3.2.4 space charge-limited conduction (SCLC)

The main objective of models based on the Space Charge-Limited conduction (SCLC) is to determine the external current of electronic carriers in a medium without traps, with only one level of trapping, or with an exponential distribution of trapping levels. The carrier concentration should be sufficiently high to produce a significant variation in the electric field and, therefore, a variation in the carriers drift velocity [START_REF] Rose | Space-Charge-Limited Currents in Solids[END_REF], [START_REF] Mark | Space-Charge-Limited Currents in Organic Crystals[END_REF].

The SCLC mechanism is similar to the transport conduction of electrons in a vacuum diode. A vacuum diode's cathode can emit electrons with a Maxwellian distribution of initial velocities (𝑉).

The corresponding charge distribution can be written by the Poisson's equation:

𝜕 2 𝑉(𝑥) 𝜕𝑥 2 = - 𝜌(𝑥) 𝜀 0 (7)
Moreover, in the steady state, the continuity equation is

𝐽 = 𝑞𝑛(𝑥)𝑣(𝑥) (8)
where

𝑣 = [ 2𝑞𝑉(𝑥) 𝑚 ] 1/2
Chapter I

35

In solid dielectrics, the SCLC is caused by the injection of electrons or holes at an ohmic contact.

The continuity equation includes the diffusion component and can be written as

𝐽 = 𝑞𝑛(𝑥)𝑣(𝑥) + 𝑞𝐷 𝑑𝑛 𝑑𝑥 (9)

I.3.2.5 Diffusion

Diffusion is a natural motion of species atoms, molecules, ions [START_REF] Müller-Plathe | No title found[END_REF] or electronic carriers [START_REF] Torpey | Double-carrier injection and recombination in insulators, including diffusion effects[END_REF] which is due to a gradient of chemical potential. At a first approximation [START_REF] Müller-Plathe | No title found[END_REF], transport by diffusion can be seen as the first derivative of species concentration in the three directions of space.

Besides, the mobility of a highly diluted ion solution is related to the diffusion coefficient by the Nernst-Einstein equation:

𝜇 = 𝑞𝐷 𝑘𝑇 (10)
This relationship presumes that ions mobility is linked to the molecular motion within the polymer [START_REF] Bamford | Ionic conductivity, glass transition, and local free volume in poly(ethylene oxide) electrolytes: Single and mixed ion conductors[END_REF]. In polymers, diffusion is often supposed negligible in the transport equations. The proposition of neglecting diffusion comes from the fact that the concentration gradients are not very significant. In fact, even when the simplest case of neutral species is considered, the diffusion coefficients can change in a significant manner, according to the size of the molecules that diffuse.

• To summarize what has been stated so far, the preceding section described numerous charge transport processes that occur in dielectrics. In our model, the hopping mechanism is the main phenomena considered in our model to explain dielectric transportation phenomena. The other mechanisms are completely ignored.

I.4 Mathematical modeling of charge transport

The bipolar charge transport (BCT) model is illustrated in Figure 12. Consider a flat sample of insulating material of thickness D that is sandwiched between a semiconducting and a stainlesssteel (i.e., anode and cathode). A positive or negative DC voltage V0 is applied to the electrode at zero time, while the other electrode is grounded. Since the electrodes radii are much greater than the dielectric sample's thickness, the edge effect can be neglected [5]. The study of charge transport in the flat sample can, therefore, be reduced to a one-dimensional domain. In such a case, most of the parameters are functions of coordinate 𝑥 along with insulation thickness and time t.

Two kinds of carriers are considered in the model, being either trapped or mobile. Charge carriers transported between shallow traps are referred to as mobile electrons and holes, whereas those being captured in the deep trapping centers are referred to as trapped electrons and holes. Next, we present the different physical phenomena that have been used in fluid modeling.

Whatever the physical model of transport is considered, and neglecting the polarization, the mathematical problem is reduced to three main equations: (11)

Where 𝑤 𝑚𝑜 𝑒,ℎ is the depth for a single trapping level, 𝑣 = 𝑘 B 𝑇 ℎ = 6.21012 × 10 -12 𝑠 -1 is the frequency of the phonons at room temperature, i.e., it is the number of times per second that the trapped electron or hole strikes the barrier of the trap and 𝑎 = 3.8 × 10 -9 𝑚 is the average distance between traps.

The motion of mobile carriers 𝑛 𝑒 and 𝑛 ℎ , results in a conductive drift current density defined as:

𝑗 𝑑𝑟𝑖𝑓𝑡 = 𝑞(𝑛 𝑒 𝜇 𝑒 + 𝑛 ℎ 𝜇 ℎ )𝐸 (15)
where 𝜇 𝑒 and 𝜇 ℎ are, respectively, electron and hole mobilities in 𝑚 2 ⋅ 𝑉 -1 . 𝑠 -1 , 𝑞 is the elementary charge in 𝐶 and 𝐸 is the applied electric field in 𝑉 ⋅ 𝑚 -1 . 𝑛 𝑒 and 𝑛 ℎ are the net charge density for electrons and holes. Hence, the transport equation is defined as:

𝑗 𝑒,ℎ (𝑥, 𝑡) = 𝑞 ⋅ 𝜇 𝑒,ℎ (𝑥, 𝑡) • 𝐸(𝑥, 𝑡) • 𝑛 𝑒,ℎ (𝑥, 𝑡) (16)

I.4.1.2 Poisson's equation

The Poisson's equation expresses the electric field distribution inside the material as a function of net charge density:

1 𝑞 • 𝜕 2 𝑉(𝑥, 𝑡) 𝜕𝑥 2 = -𝑛(𝑥, 𝑡) 𝜀 0 𝜀 𝑟 ( 17
)
Where V is the applied voltage with 𝐸 =

The source terms 𝑠 𝑒𝜇 , 𝑠 𝑒𝑡 , 𝑠 ℎ𝜇 and 𝑠 ℎ𝑡 (respectively, for mobile electrons, trapped electrons, mobile holes, and trapped holes) of Equation ( 20) elucidate and sort out the microscopic phenomena and processes inside the dielectric other than the transport process, i.e., recombination, trapping and de-trapping. The recombination terms that have been included in the model have the Langevin form [START_REF] Liu | Role of the Charge-Transfer State in Reduced Langevin Recombination in Organic Solar Cells: A Theoretical Study[END_REF], which is usually applied when at least one kind of the carriers is mobile. A reduction of the recombination rates has been considered by introducing reduction pre-factors 𝑅 1..3 .

The source term equations are defined as:

𝜕𝑛 𝑒µ 𝜕𝑡 +
The sum of the continuity equations for each type of charge leads to the global continuity equation, function of the net charge density:

𝜕𝑛(𝑥, 𝑡) 𝜕𝑡 + 1 𝑞 𝜕𝑗(𝑥, 𝑡) 𝜕𝑥 = 0 (23)
Where 𝑗 represents here the conduction current density. This equation expresses the variation in time and space of the net charge density. It is valid whatever the condition considered (transient or at equilibrium).

I.4.2 Numerical space charge models

Many numerical models to study polymeric materials performance in DC insulating systems have been proposed. The models can generally be divided into two categories, one that simulates the macroscopic behavior of the material based on conductivity data obtained as a function of temperature and electric field, and the other that describes the microscopic characteristics of space charge within a dielectric system.

The simulation results from both models will properly represent the charging density profiles and electric field distributions in the bulk of dielectrics and their progress over time compared with the experimental results. However, compared with models based on the varying conductivity, charge transport models are actually more physical since they represent the carriers transport processes in the dielectric (including charged trapping, detrapping and recombination), and especially consider the charge injection at the interface between the electrodes and the dielectric. Many researchers have further developed these kinds of models in order to achieve a better fit with experimental data, and most attempts have been made to explain the charge dynamics in solid-dielectric materials. The most famous models are represented in Table 2.

Chapter I 40 expected to recombine with their opposite species (electrons with holes). Due to the fact that oppositely charged species are being considered in the numerical computation, trapping and recombination between these species have also been taken into account in the numerical calculation.

Fukuma et al. (1994):

The model of Fukuma et al. [START_REF] Fukuma | Computer analysis on transient space charge distribution in polymer[END_REF], [START_REF] Fukuma | Numerical analysis on transient space charge distribution in XLPE[END_REF] was more comprehensive than what was previously described. The generation of charges is made by two-electrode and based on Schottky injection law. They also considered that the charges were present in the dielectric before the voltage was applied. Transport is characterized by hopping mechanism, as the two moving species can only be trapped at a low trapping level for a short time. At both interfaces, a possible barrier to carrier extraction is considered.

The recombination of mobile carriers is also considered. The barriers are symmetrical for both species (hopping, injection, and extraction).

As for the previous model, the results are space charge profiles which are compared to the experimental data of Li and Takada [START_REF] Li | Influence on spatial charge distribution of cross-linking agent residues in XLPE[END_REF]. Besides, the model has been applied for a material containing an internal interface between two layers of XLPE [START_REF] Fukuma | Numerical analysis on transient space charge distribution in XLPE[END_REF].

Chapter The model is based on the theories of Alison and Hill [START_REF] Alison | A model for bipolar charge transport, trapping and recombination in degassed crosslinked polyethene[END_REF], however detrapping is introduced here.

The approach is based on the existence of two very distinct zones in the trap energy distribution in polyethylene. The first one, which is very similar to the conduction band, has a trap depth of no more than 0,3 eV. The second varies from 0.5 to 1.5 eV due to chemical disorder.

It appears that the model is able to reproduce the essential characteristics observed under a DC stress: the space charge profiles and their dynamics under polarization and depolarization, the charge and discharge currents and electroluminescence [START_REF] Roy | Charge transport and dissipative processes in insulating polymers: experiments and model[END_REF]. this method is consistent with that obtained using a resolution in in dynamic conditions at long time, thereby validating the problem's statistical and computational treatments.

F
• The steady-state charge transport model was used to compare symmetrical and unsymmetrical transport parameters.

• The recombination terms that have been included in the model have the Langevin form.

• Highlight the importance of the role of recombination processes in the shape of the J -V characteristic. • The conventional Richardson-Schottky equation has been replaced by this injection rule since the latter causes inconsistencies in the barrier height value [START_REF] Roy | Description of bipolar charge transport in polyethylene using a fluid model with a constant mobility: model prediction[END_REF], [START_REF]Sixth International Conference on Dielectric Materials, Measurements, and Applications[END_REF]:

G. Chen et al (2009): George Chen

𝑗 𝑒,ℎ (𝑥, 𝑡) = 𝛼 exp(𝛽𝐸(𝑥, 𝑡)) (24)
Thus, 𝑗 𝑒 (0, 𝑡 ) and 𝑗 ℎ (𝐷, 𝑡 ) are the injection fluxes of electrons and holes, respectively; 𝛼 and 𝛽 are constants.

Chapter I

44

• The space charge is instantaneously defined in the space charge volume, which stays constant over time regardless of the applied field, i.e., charge transport is ignored.

• Injected or transported charges will either be trapped or recombine with the trapped charges of opposite polarity.

S. Le Roy et al (2016):

The bipolar charge transport model was extended to a cylindrical structure (cable geometry) in 2016 by Le Roy et al. [START_REF] Roy | Modelling space charge in a cable geometry[END_REF], taking into account a steady temperature gradient through the insulation. The simulation results demonstrate how geometry and temperature influence charge densities and the distributions of electrical field respectively.

A modified Schottky law was introduced at each electrode, which ensures that the injected current density for a zero field is also zero.

Y. Zhan et al (2019):

In 2019, Yupeng Zhan [START_REF] Zhan | Modelling Space Charge in HVDC Cable Insulation[END_REF] suggested a more realistic condition, i.e., thermal transients, in order to model charging behavior in a cable geometry. Based on the time dependent heat transfer equation, the temperature distribution across the insulation at each time step could be obtained. With varying temperatures and field distributions, the evolution of space charge and electric field has been studied. Furthermore, the effects of the thermal transient on the dynamics of space charging were considered by applying a current flow in cable heart. The model also includes an electric threshold field to redefine the electrode injection charges. Such that, the Schottky injection law dominates when the applied electric field exceeds the threshold, and ohmic conduction takes over under the low field.

Y. Zhan et al (2020):

A modified bipolar charge transport model is introduced by Yunpeng Zhan in 2020 [START_REF] Zhan | Space Charge Measurement and Modelling in Cross-Linked Polyethylene[END_REF], and utilized to simulate the space charge behavior in LDPE and XLPE, by considering the ionization. The used model takes into account impurity dissociation, and was able to effectively predict the features of hetero charge formation.

• The process of charge transport was defined under the constant DC field at room temperature by a constant effective mobility of carriers

• Diffusion processes have been considered with a diffusion coefficient, which is supposed to describe the natural motion of chemical species based on the space gradient of particle concentration.

Chapter I

• In this model, the charge generation in the clean LDPE is considered from injection at the electrodes following the Schottky law, and a threshold electric field has been considered (10kV/mm), at which the charge injection takes place, was introduced by [START_REF] Zhan | Modelling Space Charge in HVDC Cable Insulation[END_REF] E. Doedens et al (2020): Espen Doedens and Serdyuk [START_REF] Doedens | Space Charge Accumulation at Material Interfaces in HVDC Cable Insulation Part II-Simulations of Charge Transport[END_REF] developed a more complete model than the previously mentioned ones. This model is designed to efficiently characterize the bipolar charge transportation and space charge phenomena in solid dielectrics, particularly in XLPE and LDPE, with high DC stresses.

The following are some of the modifications made by Doedens:

• Diffusion phenomenon is considered in this model.

• Both the modified Schottky law and Fowler-Nordheim mechanism are used to describe the generation of charges on both electrodes.

• The charge injection equations have been modified to account for rough interfacial geometry (surface roughness).

• A new approach is used to reduce the injection rate even further below a particular (ohmic) field threshold. Field threshold is achieved by using a smooth Heaviside (step) function [START_REF] Doedens | Space Charge Accumulation at Material Interfaces in HVDC Cable Insulation Part II-Simulations of Charge Transport[END_REF].

• This model considers three types of traps: shallow traps, inter-level traps, and deep traps.

Which, in comparison to the previous ones, offers a more realistic situation.

• Trapping and detrapping rates are characterized by new equations considering the influence of electric field variation.

I.5 Problem formulation

One of the challenges in studying fluid models is that physical models must be very complex in order to simulate a continuously increasing number of experimental measurements, leading to an increased number of parameters, most of them are unknown in alternative ways. This type of model requires a set of initial experimental conditions, such as temperature, applied voltage, and dielectric thickness, as well as a large number of unknown parameters, such as injection barrier, mobility, trapping, and de-trapping coefficients. Most of these parameters cannot be predicted, observed, or estimated by independent experiments. In other words, it is an elusive task to compute this kind of parameters that best fit the experimental data. A variety of methods have been used to approximate the parameter values, some of which are based on ab initio methods based on the materials used, while others are based on the charge density results obtained by PEA. Other researchers used an optimization algorithm to approximate these parameters by solving the inverse problem using PEA's experimental data. The following are several parameter approximations made by various researchers:

• The PEA measurements have been used by several researchers [1], [START_REF] Perlman | Drift mobility determinations using surface-potential decay in insulators[END_REF], [START_REF] Mazzanti | A space-charge based method for the estimation of apparent mobility and trap depth as markers for insulation degradation-theoretical basis and experimental validation[END_REF], to approximate the charge mobility. Depending on the type of carrier and the temperature, the final parametrization produces mobilities in the range of 10 -12 -10 -14 𝑚 2 𝑉 -1 𝑠 -1 .

• Quirke et al. [START_REF] Meunier | Molecular modeling of electron traps in polymer insulators: Chemical defects and impurities[END_REF] proposed a trapping depth range that is close to 1 eV using ab initio methods. Other research on LDPE and XLPE [START_REF] Mazzanti | A space-charge based method for the estimation of apparent mobility and trap depth as markers for insulation degradation-theoretical basis and experimental validation[END_REF], in particular, has found the similar order of trapping depth parameters (0.9-1 eV), which correspond to chemical traps (double bonds).

• For electron injection, the theoretical injection barrier height corresponding to a goldpolyethylene interface is 5 eV. When this value is applied to the Schottky law, the simulated current density is ten decades lower than the experimental value. This theoretical value does not account the local interface conditions (local field strengthening on specific locations, chemical impurities forming deep traps) and therefore cannot be used in simulations. Other research proposed different values of the injection barriers of the order of 1 eV based on the outcomes of charge transport models in polyethylene [START_REF] Fukuma | Computer analysis on transient space charge distribution in polymer[END_REF],

[51], [START_REF] Hozumi | Space charge behavior in XLPE cable insulation under 0.2-1.2 MV/cm dc fields[END_REF].

• Some of the parameters, such as recombination coefficients, could not be defined easily by experiments. For this purpose, the criticality of the parameters was examined, i.e., by studying the effects of modifying the given parameter on the model outputs [START_REF] Roy | Numerical model for studying dynamic space charge behavior in polyethylene[END_REF]. It has been shown [START_REF] Roy | Numerical model for studying dynamic space charge behavior in polyethylene[END_REF] that the recombination coefficients do not play a major role in simulating results at low temperatures or short periods of polarization. By considering one of these cases, recombination coefficients do not change the current density to a significant degree, and they weakly affect the charge density in the dielectric, in the limit of the investigated range of recombination rates from 10 -5 to 1 m 3 /Cs.

To summarize what has been presented so far, many alternative methods and techniques have been utilized to find an adequate set of parameters that can match any experimental data obtained by PEA or current measurements. However, the optimal values of the parameters have not yet been achieved.

Chapter I

47

As a conclusion, the major objective of this study is to provide an original approach that contributes to identifying a unique set of parameters that provide a good correlation between experimental and simulated results of charge and current densities using any experimental conditions.

II.1 Introduction

Optimization is a branch of applied mathematics concerned with minimizing a function in a defined domain while keeping various variables (or unknowns) in consideration [1]- [3]. In other words, optimization aims to identify values for the variables that will best achieve the optimum value of the function. The function is known as the cost function or the objective function in this context. The objective or cost function is a mathematical equation describing a specific problem that should be minimized or maximized. For example, the relative error between experimental data and a numerical model is an example of a cost function that should be minimized to fit the model to the experimental measurements.

Historically, optimization approaches were initially used to solve difficulties with human logistics and transportation management [4]- [6]. The discipline of optimization is still a very active area of research. For example, optimization tools are used in engineering and manufacturing to find the best and most durable material structures [7]. Besides, while speaking about business, optimization is usually used in business domain to find the highest profits with the lowest costs [8].

In recent years, the field of optimization has gotten much interest, owing to significant advances in computer technology, such as the availability of efficient software, high-speed computers, and artificial neural networks. Optimization is present in every aspect of life, such as, airline scheduling, finance, internet routing and engineering design [9], [10]. We continuously attempt to optimize anything in engineering and industry, reducing costs and energy consumption or increasing profit, production, performance, time, and efficiency.

The objective function cannot be described in clear analytical terms (i.e., this sort of functions is called black-box functions) in most engineering and industrial applications since the dependency on the design factors of the objective are complex and implicit. Therefore, this form of black-box optimization frequently necessitates the use of numerical modeling, which is often computationally costly, such as computational fluid dynamics [11] or finite element analysis [12].

In addition, all optimization strategies are iterative, requiring several function evaluations.

Accordingly, every method that enhances simulator efficiency or minimizes the number of function evaluations is essential. Indeed, in numerical optimization, it is necessary to properly characterize the problem being addressed to determine which algorithm to employ.

Chapter II

56

There is no universal optimization technique; instead, various algorithms are customized to a particular type of optimization problem. It is typically the user's responsibility to choose the best algorithm for a specific application. This is a critical decision since it might determine whether the problem is addressed quickly or slowly and whether the optimal solution can be found or not. We must choose an optimization algorithm that can effectively find a solution once applied to the model. Good mathematical formulas known as optimality criteria may be employed in numerous situations to verify that the current set of variables reflects the problem's solution. If the optimality or stopping criteria aren't achieved, optimality criteria might provide valuable contribution into enhancing the existing estimated solution. Techniques like sensitivity analysis (e.g., Sobol Index) [13], which shows the inputs sensitivity on the outputs of a model, can improve the model's

behavior by fixing or eliminating some parameters that do not influence the model outputs.

This chapter aims to offer a detailed discussion of numerical optimization methods and a decision tree of optimization approaches that the reader may use to understand each technique employed throughout the chapter.

II.2 Optimization basic concepts

Optimization is the mathematical term for the minimization or maximizing of a function subject to variable constraints. The optimization problem is stated as follows:

min 𝑥∈ℝ 𝑛 𝑓(𝑥) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑥 ∈ 𝒳 (1)
Where 𝑓: ℝ 𝑛 → ℝ is the objective function, which is a (scalar) function in terms of 𝑥 that we aim to maximize or minimize; 𝑥 is the vector of 𝑛 variables, which are also known as unknowns or parameters:

𝑥 = [𝑥 1 , 𝑥 2 , … , 𝑥 𝑛 ]
The elements in this vector can be modified in order to minimize the objective function 𝑓 . A solution or minimizer is denoted by 𝑥 ⋆ from among all points in the feasible set 𝒳 that minimizes the objective function, such that, 𝑓(𝑥 ⋆ ) ≤ 𝑓(𝑥), ∀𝑥 ∈ 𝒳.

II.2.1 Minimum and maximum

The maxima and minima (the corresponding plurals of maximum and minimum) of a function are

Since the search for a maximum is equivalent to the search for a minimum, we will limit ourselves in this chapter to minimization problems. Furthermore, if 𝒳 = ℝ 𝑛 , we claim the issue is unconstrained. If 𝒳 is a closed set of ℝ 𝑛 , then we say that the problem is constrained.

II.2.2 Gradient and Hessian

The gradient and the hessian of 𝑓 are the vectors and matrices defined by:

𝛻𝑓(𝑥) = ( 𝜕𝑓(𝑥) 𝜕𝑥 𝑖 ) 𝑖=1..𝑛 = ( 𝜕𝑓(𝑥) 𝜕𝑥 1 𝜕𝑓(𝑥) 𝜕𝑥 2 ⋮ 𝜕𝑓(𝑥) 𝜕𝑥 𝑛 ) ∈ ℝ 𝑛 (3) 𝛻 2 𝑓(𝑥) = ( 𝜕 2 𝑓(𝑥) 𝜕𝑥 𝑖 𝜕𝑥 𝑗 ) 𝑖,𝑗=1..𝑛 = ( 𝜕𝑓 2 (𝑥) 𝜕𝑥 1 2 𝜕𝑓 2 (𝑥) 𝜕𝑥 1 𝜕𝑥 2 ⋯ 𝜕𝑓 2 (𝑥) 𝜕𝑥 1 𝜕𝑥 𝑛 𝜕𝑓 2 (𝑥) 𝜕𝑥 2 𝜕𝑥 1 𝜕𝑓 2 (𝑥) 𝜕𝑥 2 2 ⋯ 𝜕𝑓 2 (𝑥) 𝜕𝑥 2 𝜕𝑥 𝑛 ⋮ ⋮ ⋱ ⋮ 𝜕𝑓 2 (𝑥) 𝜕𝑥 𝑛 𝜕𝑥 1 𝜕𝑓 2 (𝑥) 𝜕𝑥 𝑛 𝜕𝑥 2 ⋯ 𝜕𝑓 2 (𝑥) 𝜕𝑥 𝑛 2 ) ∈ ℝ 𝑛×𝑛 (4)
For a matrix 𝐴 ∈ ℝ 𝑛×𝑛 , we recall the definitions of positivity (noted ≻ 0) and semi-definite positivity (noted ≽ 0):

We say that 𝐴 ≽ 0 if and only if:

𝑣 𝑇 𝐴𝑣 ≥ 0, ∀𝑣 ∈ ℝ 𝑛 (5)
And 𝐴 ≻ 0 if and only if:

𝑣 𝑇 𝐴𝑣 > 0, ∀𝑣 ∈ ℝ 𝑛 \{0 ℝ 𝑛 } (6)

II.2.3 Critical Points

Figure 14 shows a univariate function 𝑓(𝑥) with several critical points where the derivative is zero, which are helpful in explaining optimization concerns. When minimizing 𝑓(𝑥), we would like to find a global minimizer or a value of 𝑥 for which 𝑓(𝑥) is minimized. Unfortunately, proving that a specific point is at a global minimum has often been challenging. Usually, the most we can do is to check whether it is at a local minimum or not. The different types of critical points are defined as:

• Local: A point 𝑥 ⋆ is a weak local minimizer if: ∃δ > 0 such that 𝑓(𝑥 ⋆ ) ≤ 𝑓(𝑥), ∀ 𝑥 ⋆ ∈ 𝐵(𝑥, δ).
Where

𝐵(𝑥, δ) = {𝑥 ∈ ℝ 𝑛 s. t. ‖𝑥 -𝑥 ⋆ ‖ < δ }
A strong local minimizer, also known as a strict local minimizer, is a point that strictly minimizes 𝑓 within a neighborhood 𝐵(𝑥, δ) (i.e., a ball of center 𝑥 and radius δ). In other words, 𝑥 ⋆ is a strong local minimizer if:

∃δ > 0 such that 𝑓(𝑥 ⋆ ) < 𝑓(𝑥), ∀ 𝑥 ⋆ ∈ 𝐵(𝑥, δ).
A weak local minimizer is a local minimizer that is not a strong local minimizer. Figure

II.2.4 Necessary Conditions for Local Minima

Many approaches of numerical optimization target local minimum requirements. We usually don't know if the local minimum is a global minimum or not. The conditions we discuss in this section assume that the objective function is differentiable. We also assume that the problem is unconstrained in this part.

The following conditions are necessary for 𝑥 to be at a local minimum of 𝑓:

• ∇𝑓(𝑥) = 0, the first-order necessary condition (FONC)

• ∇ 2 𝑓(𝑥) is positive semi-definite, the second-order necessary condition (SONC)

A simple analysis may be used to determine the FONC and SONC. 𝑥 ⋆ must be less than the values surrounding it in order to be at a local minimum:

𝑓(𝑥 ⋆ ) ≤ 𝑓(𝑥 + ℎ) ⟺ 𝑓(𝑥 ⋆ ) -𝑓(𝑥 + ℎ) ≤ 0 (7)
If we write the second-order approximation for 𝑓(𝑥 ⋆ ), we get:

𝑓(𝑥 ⋆ + ℎ) = 𝑓(𝑥 ⋆ ) + ∇𝑓(𝑥 ⋆ ) ⊤ ℎ + 1 2 ℎ ⊤ ∇ 2 𝑓(𝑥 ⋆ )ℎ + O(ℎ 3 ) (8)
We know that at a minimum, the first derivative must be zero, and we neglect the higher-order terms. Rearranging, we get:

1 2 ℎ ⊤ ∇ 2 𝑓(𝑥 ⋆ )ℎ = 𝑓(𝑥 + ℎ) -𝑓(𝑥 ⋆ ) ≥ 0 (9)
This is the definition of a positive semi-definite matrix, and we recover the SONC. The FONC and SONC are essential for optimality, but they are not sufficient. For unconstrained optimization of a twice-differentiable function, if the FONC is met and ∇ 2 𝑓(𝑥) is positive definite, 𝑥 is guaranteed to be at a strong local minimum. These conditions are collectively known as the second-order sufficient condition. Figure 15 shows some examples about SONC and FONC of the necessary but insufficient conditions for strong local minima.

II.3 Types of Optimization Problems

The general form problems (1) may be classified in terms of the nature of the objective function and the nature of the constraints. Chapter II 61

II.3.1 Constrained and unconstrained optimization problems II.3.1.1 Unconstrained optimization problems

As previously mentioned, unconstrained optimization problems arise when 𝒳 = ℝ 𝑛 in (1). That is to say, unconstrained optimization are the problems when we try to minimize an objective function that depends on real variables with no constraints on their values. An unconstrained optimization problem can be of any kind (linear or nonlinear). The mathematical formulation is defined as:

min 𝑥 𝑓(𝑥), (10)
Where 𝑥 ∈ ℝ n is a real vector of n ≥ 1 components and 𝑓: ℝ n → ℝ is the cost function.

Even with some natural variables limitations, it can be safe to ignore the constrained as they do

not influence the answer and do not affect the algorithm's behavior.

In order to solve problem 12, suppose that we know how to evaluate (directly or by finite difference) the gradient ∇𝑓(𝑥) of 𝑓 at any point of the search space. Descent methods are iterative algorithms which starts from initial point 𝑥 0 , construct a sequence of iterates 𝑥 𝑘 defined as:

𝑥 𝑘+1 = 𝑥 𝑘 + 𝑠 𝑘 𝑑 𝑘 ,
Where 𝑠 𝑘 is the step length of the displacement, 𝑑 𝑘 ∈ ℝ n is the descent direction of 𝑓 at 𝑥 𝑘 verifying the inequality:

𝛻𝑓(𝑥) 𝑑 𝑘 < 0,
A descent algorithm is thus completely characterized by the way the descent directions 𝑑 𝑘 are constructed and the step length 𝑠 𝑘 that one takes in these directions. Furthermore, the descent .

II.3.1.2 Constrained optimization problems

Constrained optimization problems [14]- [16] are addressed in various scientific domains, such as structural optimization, engineering design, economics, allocation, and location concerns.

Constrained optimization (CO) problems arise from models in which constraints play an essential role, such as imposing financial limitations in an economic problem or shape constraints in a design problem.

The CO problem can be represented as the following nonlinear programming problem:

min 𝑓(𝑥) , 𝑥 ∈ 𝒳 ⊂ ℝ 𝑛 , (11)
Subject to the linear or nonlinear constraints:

𝑔 𝑖 (𝑥) = 0, 𝑖 = 1, … , 𝑚. (12)
ℎ 𝑖 (𝑥) ≤ 0, 𝑖 = 1, … , 𝑚.

Constraints can also be expressed as range constraints, which specify that the values of some decision variables must lie within specific closed intervals of ℝ, such as 𝑎 ≤ 𝑥 ≤ 𝑏 with 𝑎, 𝑏 ∈ ℝ .Besides, a and b can be unbounded such that [𝑎, 𝑏] = ] -∞, 𝑏] or [𝑎, +∞ [ (i.e., negative or positive constraints). Another well-known type of constraints known as nonlinear inequalities represents the complex connections between the variables. For example, the polynomial constraints such as ‖𝑥‖ 2 2 ≤ 𝑎.

The CO problems can be solved using several methods, such as: Lagrange multipliers, Penalty methods, etc...

• Lagrange multipliers:

The method of Lagrange multipliers [17], [18] is a strategy for finding the local maxima and minima of a function subject to equality constraints. The main concept is to transform a constrained problem into a form such that the derivative test of an unconstrained problem can be applied. The connection between the function's gradient and the gradients of the constraints easily leads to a reformulation of the original problem known as the Lagrangian function.

The approach is summarized as follows: in order to get the maximum or minimum of a function 𝑓(𝑥) subjected to the equality constraint 𝑔 𝑖 (𝑥), form the Lagrangian function,

Chapter II 63 ℒ(𝑥, 𝜆) = 𝑓(𝑥) -𝜆𝑔(𝑥) (14) then calculate the stationary points of ℒ considered as a function of 𝑥 and the Lagrange multiplier 𝜆. The solution to the initial constrained optimization is always a saddle point of the Lagrangian function, which may be found among the stationary points by the definiteness of the bordered Hessian matrix. The Lagrange multiplier approach is widely applied to solve complicated constrained optimization problems. Furthermore, the Lagrange multiplier approach is generalized by the Karush-Kuhn-Tucker conditions, which may account for inequality constraints of the form ℎ(𝑥) ≤ 0.

• Penalty methods: A penalty approach [19] replaces a constrained optimization problem with a set of unconstrained problems whose solutions converge to the solution of the original constrained problem. Unconstrained problems are obtained by adding a term to the objective function known as penalty function, which consists of a penalty parameter multiplied by a measure of constraint violation. When the constraints are violated, the measure of violation is nonzero, and it is zero in the region where the constraints are not violated.

The method can be summarized as follows: find the maximum or minimum of a function 𝑓(𝑥) subjected to the equality constraint ℎ 𝑖 (𝑥) ≤ 0.

This problem can be solved as a series of unconstrained minimization problems:

min 𝐶 𝑘 (𝑥) = 𝑓(𝑥) + 𝜎 𝑘 ‖𝐼(ℎ 𝑖 (𝑥))‖ 1 2 , ( 15
)
where 𝐼(ℎ 𝑖 (𝑥)) = max(0, ℎ 𝑖 (𝑥)) 2 .

In the above equations, 𝐼(ℎ 𝑖 (𝑥)) is the exterior penalty function, while 𝜎 𝑘 are the penalty coefficients. In each iteration 𝑘 of the method, we increase the penalty coefficient 𝜎 𝑘 . Solve the unconstrained problem (17) and use the solution as the initial guess for the next iteration. Solutions of the successive unconstrained problems will eventually converge to the solution of the original constrained problem.

II.3.2 Linear and nonlinear programming problem II.3.2.1 Linear programming problems

The problem is called a linear programming problem [3], [15] when the objective function and all constraints are all linear functions in terms of 𝑥. These problems are the most commonly

Chapter II 64 formulated and solved optimization problems with management, finance, and economics applications.

A problem is considered to be linear when the objective and constraints are linear. Every linear problem may therefore be placed as follows:

min 𝑐 𝑇 𝑥 s.t { 𝐴𝑥 = 𝑏 𝑥 ≥ 0 (16)
Where 𝑐 ∈ ℝ n , 𝐴 ∈ ℝ p×n , 𝑎𝑛𝑑 𝑏 ∈ ℝ p .

The reason for the popularity of linear forms for constraints and objectives in problem formulation is that they are the simplest algorithms to be used.

II.3.2.2 Nonlinear programming problems

Nonlinear programming problems (NLP) [3], [20], in which at least part of the constraints or the objective are nonlinear functions, tend to arise naturally in the physical sciences and engineering, and are becoming frequent in management and economic sciences as well. Nonlinear programming is certainly a helpful method, since numerous aspects of our world do not behave linearly.

In NLP, a nonlinear objective function may be minimized or maximized subject to bounded constraints, linear constraints, or nonlinear constraints, where the constraints might be inequalities or equalities. A nonlinear minimization problem is an optimization problem of the form:

min 𝑓(𝑥) 𝑠. 𝑡. 𝑔 𝑖 (𝑥) ≤ 0 ∀𝑖 ∈ {1, … , 𝑚} ℎ 𝑗 (𝑥) = 0 ∀𝑗 ∈ {1, … , 𝑝} (17)

II.3.3 Least Squares optimization

Least Square Optimization (LSO) [1], [21] is one of the most popular optimization problems, and Chapter II

65

One of the advantages of using LSO is the ability to fit a wide variety of functions is the most significant benefit of nonlinear least squares optimization over many other approaches. While many scientific and technical processes may be successfully represented using linear models or other relatively basic types of models, many others are fundamentally nonlinear.

Theory: Consider a set of 𝑚 data points (𝑥 1 , 𝑦 1 ), (𝑥 2 , 𝑦 2 ), … , (𝑥 𝑚 , 𝑦 𝑚 ), where 𝑥 i is an independent variable and 𝑦 i is a dependent variable whose value is found by observation. On the other side, we have the model function 𝑦 =𝑓(𝑥,𝑃), that in addition to the variable 𝑥 also depends on n parameters, such that,

𝑃 = (𝑃 1 , 𝑃 2 , … , 𝑃 𝑛 )
The main target of the least-squares method is to identify a point 𝑃 ⋆ which minimizes the sum of squares of nonlinear function 𝐶(𝑃).

A model's fit to a data point is quantified by its residual, which is defined as the difference between the experimental value and the value predicted by the model:

𝑟 𝑖 (𝑃) = 𝑦 𝑖 -𝑓(𝑥, 𝑃),
The least-squares method determines the optimal parameter values by minimizing the residual's sum of squared , S:

𝐶(𝑃) = ∑ 𝑟 𝑖 2 𝑛 𝑖=1 = 1 2 𝑟 𝑇 (𝑃)𝑟(𝑃), 𝑃 ∈ ℝ 𝑛 , 𝑚 ≥ 𝑛 ∇𝐶(𝑃) = 𝐽 𝑟 (𝑃)𝑟(𝑃) ∇ 2 𝐶(𝑃) = 𝐽 𝑟 (𝑃) ⊤ 𝐽 𝑟 (𝑃) + 𝐻 𝑟 (𝑃)𝑟(𝑃)
Where 𝐽 𝑟 (𝑃) and 𝐻 𝑟 (𝑃) are respectively the Jacobian and Hessian of 𝑟.

This Hessian is usually very expensive to calculate. In addition, in the neighborhood of the optimum (theoretical validity domain of Newton's method) 𝑟 will be small, making the contribution of the term 𝐻 𝑟 (𝑃)𝑟(𝑃) negligible. It is therefore wise to ignore it and thus consider that:

∇ 2 𝐶(𝑃) = 𝐽 𝑟 (𝑃) ⊤ 𝐽 𝑟 (𝑃)
This approximation is known as the Gauss-Newton hypothesis [22]. It allows to approach ∇ 2 𝐶(𝑃)

by first derivatives and to ensure, on the other hand, that ∇ 2 𝐶(𝑃) is always positive semi-definite.

This hypothesis leads to the so-called Gauss-Newton algorithm.

II.4 Based derivative and free derivative II.4.1 Derivative-Based Algorithms

The information of derivatives is used in derivative-based or gradient-based algorithms [20]. The objective function must be sufficiently smooth to allow for the existence of first (and typically second) derivatives. For smooth continuous-domain problems, derivative-based optimization is effective in locating local optima. However, it may cause some issues when 𝑓 is disconnected or when 𝑓 is costly to evaluate, non-smooth, or noisy, causing derivatives to be useless. Derivativebased algorithms are frequently employed in various applications and discrete modeling [23], [24].

The problem to be solved is to optimize an objective function 𝑓:

min 𝑥∈ℝ 𝑛 𝑓(𝑥) 𝑠𝑢𝑏𝑗𝑒𝑐𝑡 𝑡𝑜 𝑥 ∈ 𝒳 (18)
The variables are updated at each major iteration 𝑘 using

𝑥 𝑘+1 = 𝑥 𝑘 + 𝛼 𝑘 𝑑 𝑘
Where 𝑑 𝑘 is the search direction for major iteration 𝑘, and 𝛼 𝑘 is the accepted step length from the line search.

Thus, to solve problem (18), one must solve the following two subproblems:

1) The computation of a search direction 𝑑 𝑘 , Algorithm 2: Gauss-Newton

Input condition: 𝑃 0 ∈ ℝ 𝑛 2) Then, search for an acceptable step size 𝛼 𝑘 .

Derivative-based algorithms use derivative information of 𝑓 to find a good search direction, since for example the gradient gives the direction of steepest descent. The various types of derivativebased algorithms are classified based on the method that is used for computing the search direction.

Some difficulties that may affect the derivative-based algorithms convergence are:

• No analytical description of the function (e.g., simulation).

• Multiple global optima (e.g., multimodal).

• Stochastic function evaluation (e.g., noisy).

• Discontinuous objective function (e.g., regions with invalid solutions).

Two of the most important algorithms are given as an example of derivative based algorithms:

Gradient descent and Newton's methods.

II.4.1.1 Gradient descent methods (First-order algorithms)

Gradient methods are based on the following idea: since the gradient 𝛻𝑓(𝑥 𝑘 ) is the direction of the largest increase of 𝑓 at the point 𝑥 𝑘 , it is natural to move in the opposite direction. Thus, the directions of descent are defined by 𝑑 𝑘 = -𝛻𝑓(𝑥 𝑘 ) , generating a sequence of iterates defined by:

𝑥 𝑘+1 = 𝑥 𝑘 + 𝛼 𝑘 𝛻𝑓(𝑥 𝑘 )
‖𝛻𝑓(𝑥 𝑘 )‖ , with 𝑠 𝑘 > 0

The distinction between all approaches based on gradient descent methods will be in the choice of the step length 𝛼 𝑘 .

Algorithm 3: Gradient descent

Input condition: 𝑥 0 ∈ ℝ 𝑛 The general principle of this method is to minimize about 𝑥 𝑘 , the quadratic approximation of 𝑓, 𝑚 𝑘 :

𝑚 𝑘 (𝑥 𝑘 + 𝑠) = 𝑓(𝑥 𝑘 ) + ∇𝑓(𝑥 𝑘 ) ⊤ ⋅ 𝑠 + 1 2 𝑠 ⊤ ⋅ 𝛻 2 𝑓(𝑥 𝑘 ) ⊤ ⋅ 𝑠
The minimum 𝑠 of 𝑚 𝑘 , called Newton's direction at point 𝑥 𝑘 can be found by solving:

∇𝑚 𝑘 (𝑥 𝑘 + 𝑠) = 0
Thus,

𝛻 2 𝑓(𝑥 𝑘 ) ⋅ 𝑠 𝑘 = -∇𝑓(𝑥 𝑘 )
An iteration of Newton's algorithm is then written:

𝑥 𝑘+1 = 𝑥 𝑘 + 𝑠 𝑘

II.4.2 Derivative free Algorithms

Derivative-free optimization [25] is a branch of mathematical optimization that does not require derivative information to identify optimal solutions. Sometimes information regarding the derivative of the objective function 𝑓 is unavailable, unreliable, or impractical to obtain. For example, 𝑓 might be noisy, non-smooth, or time-consuming to evaluate, so that methods that rely It is well known that the derivatives of any function one wants to optimize, provide a lot of important information. That is, for continuously differentiable functions, the "standard" mathematical characterization of a local minimum, provided by the first-order necessary conditions, necessitates that the first-order derivatives be zero. However, there have always been many instances when derivatives are unavailable or inaccurate for several reasons.

Indeed, optimization with no derivatives is viewed as one of the most significant and challenging fields and one with huge practical promise in computer science and engineering. Derivative-free optimization is now in high demand due to the increasing complexity of mathematical modeling, the sophistication of scientific computing.

In many optimization problems arising from engineering, scientific, and artificial intelligence applications, objective functions, and constraints are available only as the output of a black-box or simulation model that does not provide derivative information. Such situations demand the use of derivative-free methods.

Derivative-based algorithms are efficient, but they may impose certain strict requirements on the objective functions. If the objective functions are discontinuous, derivative-free methods may be more efficient and logical. One of the most prevalent reasons for employing derivative-free approaches in the early days of nonlinear optimization was probably the user's lack of knowledge.

Users understood they wanted to enhance their present solution, but they wanted to use something simple and easy to explain, so they utilized free derivative approaches, even when better algorithms were available.

II.5 Conclusion

Numerical optimization is not widely utilized in our community of research. For this reason, this chapter reviewed the basic concepts of optimization that should be well understood in order to know how to apply optimization on a specific problem. The most important types of optimization problems depending on the type of the cost functions and constraints were presented and discussed.

All the information provided in this chapter will be used in the following chapter to identify the nature of our problem, which will help in selecting the most suitable algorithms that could help in solving our problem.

III.1 Introduction

Fitting a model to a set of experimental observations is an important topic that is addressed by various researchers in many scientific fields. Consider a set of experimental data {(𝑥 1 , 𝑦 1 ), (𝑥 2 , 𝑦 2 ), … , (𝑥 𝑛 , 𝑦 𝑛 )} ∈ ℝ 2𝑛 , which relate an input variable 𝑥 with an output variable 𝑦.

The unknown relationship between 𝑥 and 𝑦 can be modeled by a functional relationship 𝑦 = 𝑓(𝑥, 𝑃), where 𝑃 are unknown parameters related to the model and 𝑦 is the output produced by the model. Thus, the problem of modeling the unknown relationship between 𝑥 and 𝑦 can be reduced to finding a value of 𝑃 that minimizes the error measure between model outputs 𝑓(𝑥 𝑖 , 𝑃)

and experimental data 𝑦 𝑖 . In the Nonlinear Least Squares (NLS) problem, the error measure is considered as the sum of the squares of the deviations between experimental data and simulation data. Thus, the methodology is as follows: assume that we need to find a set of parameters vector 𝑃 ⋆ such that it minimizes the cost function 𝐶(𝑃), defined as:

𝐶(𝑃) = ‖𝑦 -𝑓(𝑥, 𝑃)‖ 2 2 (19)
The Nonlinear Least Squares (NLS) problem to be solved is thus defined as:

min 𝑃∈ℝ 𝑛 𝐶(𝑃) 𝑠. 𝑡. 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏 (20)
Where 𝑙𝑏 and 𝑢𝑏 represents the upper and lower bounds of the parameters, respectively;

In our case, the bipolar charge transport model is considered to be optimized. Such models require some Experimental Conditions (EC) such as temperature, applied voltage, polarization time, etc., alongside a set of parameters such as injection barriers, mobilities, trapping, de-trapping, and recombination coefficients. Most of these parameters are not predictable, observed, or estimated from independent experiments. This makes the simulation of experimental data liable to ambiguities. The main target of this research is to find a unique set of parameters that provides a good correlation between experimental (𝑓 𝑒𝑥𝑝 ) and simulated data (𝑓 𝑠𝑖𝑚 (𝑃)) for any EC (Figure 16). The outputs of the BCT model could not be described by an analytic function (e.g., simulation). Thus, the derivative cannot be calculated directly; instead, it should be approximated.

Chapter III 76 Thus, the BCT model is considered as a blackbox function. Several optimization algorithms are utilized to compute the optimal set of parameters related to the bipolar charge transport model. Iterative optimization methods start with a guess for the variable 𝑃 and iterates through a series of better guesses (called "iterations") until they reach a minimum. One algorithm differs from another based on the technique of moving from one iteration to the next. In our case, the cost function has no analytical description; hence, the first and second derivatives are obtained by an approximation instead of being analytically computed (i.e., the first and second derivatives are approximated using finite difference methods).

Effective algorithms should have the following characteristics:

• Robustness: The algorithm should find an appropriate optimal solution for any starting point ranging in the domain of search.

• Efficiency: The time of computation, number of iterations, and the functions evaluation should not be very big.

• Accuracy: They should be able to identify a solution with good accuracy.

Even if an efficient optimization technique is used, evaluating the objective functions might take a long time based on the cost function complexity. For the BCT model, some simulations might take hours, days, or even weeks. This means that finding the most suitable optimization algorithm for solving our problem is a challenging task.

This chapter provides a detailed investigation of optimality requirements, as well as accuracy and efficiency analysis that demonstrates the strong and weak points of the most significant algorithms that could be applied for optimizing the BCT model. Five optimization algorithms are considered for estimating the unknown set of parameters related to the BCT model. These algorithms are defined as: Trust Region Reflective Algorithm (TRRA), Levenberg-Marquardt (LM), Nelder-Mead (NM), Genetic Algorithm (GA), and Particle Swarm (PS). These algorithms will be employed to estimate a set of 8 unknown parameters that will be defined in the following sections.

All the considered optimization algorithms will be compared in detail by examining the robustness, efficiency, and accuracy of each one when applied to optimize the BCT model.

III.2 Algorithms

With the increasing reliance on modeling optimization problems in practical applications, several

III.2.1 Direct search Algorithms

Direct optimization algorithms are used for objective functions for which the gradient cannot be directly calculated. Gradient information is approximated from the outputs of the objective function using the finite difference methods. The algorithms have deterministic procedures and often assume that the objective function has a single global optimum.

There are many examples of direct search algorithms, however, the most suitable algorithms for solving our problem will be presented and used in this chapter. They are: Trust Region Reflective Algorithm (TRRA), Levenberg-Marquardt algorithm (LM), and Nelder-Mead algorithm (NM).

III.2.1.1 Trust Region Reflective Algorithm (Constrained)

TRRA [1] is one of the most important numerical optimization methods in solving nonlinear Least squares (NLS) problems. The TRRA is classified as one of the derivative-based algorithms where it uses both gradient and hessian information to find a good search direction.

Assume that we are interested in finding the value of a set of parameters vector of real variables 𝑃 ⋆ such that it minimizes the cost function 𝐶(𝑃):

min 𝑃∈ℝ 𝑛 𝐶(𝑃) 𝑠. 𝑡. 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏 (21)
The concept of the trust region method is to approximate 𝐶(𝑃) with a quadratic function 𝑚 𝑘 [2],

which reflects the behavior of function 𝐶(𝑃) in a neighborhood of 𝑃, defined as 𝑁(𝑃), which is called the ''trust region'', around the current point 𝑥. Following Taylor's theorem, 𝐶(𝑃 𝑘 ) can be approximated by the quadratic function 𝑚 𝑘 (𝑠), in the form:

𝑚 𝑘 (𝑃 𝑘 + 𝑠) = 𝐶(𝑃 𝑘 ) + ∇𝐶(𝑃 𝑘 ) ⊤ ⋅ 𝑠 + 1 2 𝑠 ⊤ ⋅ 𝛻 2 𝐶(𝑃 𝑘 ) ⊤ ⋅ 𝑠 (22)
where ∇𝐶(𝑃 𝑘 ) ∈ ℝ 𝑛 is the gradient of 𝐶(𝑃) and ∇ 2 𝐶(𝑃 𝑘 ) ⊂ ℝ 𝑛×𝑛 the hessian of 𝐶(𝑃) evaluated at 𝑃 𝑘 .

For instance, there are directions along which 𝑚 𝑘 (𝑠) is unbounded from below. In this case, ‖𝑠‖ is infinite. Therefore, to ensure the convergence of the iteration, the TRRA defines a region around the iterate 𝑃 𝑘 in which we can trust 𝑚 𝑘 , called "Trust Region". In practice, the trust region N is usually considered as a sphere defined by Euclidean norm ‖𝑠‖ ≤ 𝑟 𝑘 , where the scalar 𝑟 𝑘 is called the trust region radius. The two key points of trust region algorithms are therefore computing and updating the radius of the trust region 𝑟 𝑘 and solving the trust region sub-problem. The trust region sub-problem to be solved is thus defined as:

min 𝑠∈ℝ 𝑛 𝑚 𝑘 (𝑃 𝑘 + 𝑠) 𝑠. 𝑡. ‖𝑠‖ 2 ≤ 𝑟 𝑘 (23)
If 𝐶(𝑃 𝑘 + 𝑠) < 𝐶(𝑃 𝑘 ), the current point 𝑃 𝑘 is updated to be 𝑃 𝑘 + 𝑠, and the trust region 𝑁(𝑃) can remain the same for the next step. Otherwise, 𝑃 remains unchanged and trust region 𝑁(𝑃) will be reduced for the next iteration and the trust region sub-problem will be solved again. There are

Chapter III 79 several approximation methods for solving the trust region sub-problem, where the two most used are the Dogleg method [3] and the of Seihaug-Toint [4].

The radius of the trust region is determined by trial and error. At the first iteration, 𝑟 0 is set arbitrarily. Then, at each iteration, the quality of the solution 𝑠 𝑘 ⋆ of the sub-problem is evaluated and 𝑟 𝑘 is adjusted according to this evaluation. The evaluation of the quality of the solution is carried out using the quantity:

𝑙 = 𝐶(𝑃 𝑘 + 𝑠 𝑘 ⋆ ) -𝐶(𝑃 𝑘 ) 𝑚 𝑘 (𝑃 𝑘 + 𝑠 𝑘 ⋆ ) -𝑚 𝑘 (𝑃 𝑘 )
which corresponds to the ratio of the increase (positive or negative) of 𝐶(𝑃) to quadratic approximation function 𝑚 𝑘 . If 𝑚 𝑘 is reliable, this amount should be close to, or even greater than 1. Inversely if it is close to 0, 𝑚 𝑘 is not reliable. Practically, in order to characterize these cases and to update 𝑟 𝑘 , let us introduce the constants 𝜂 1 and 𝜂 2 such that 0 < 𝜂 1 ≤ 𝜂 2 < 1. Defines three cases:

If 𝑙 ≥ 𝜂 2 then 𝑚 𝑘 is very reliable and 𝑟 𝑘 is doubled

If 𝜂 1 ≤ 𝑙 < 𝜂 2 then the match between 𝑚 𝑘 and 𝐶 is not perfect, but if it allowed to reduce the value of 𝐶, then it can be assumed reliable and 𝑟 𝑘 is not modified.

If 𝑙 < 𝜂 1 then 𝑚 𝑘 is not reliable and 𝑟 𝑘 is reduced to

1 2 ‖𝑠 𝑘 ⋆ ‖.
To give an idea, we fix in practice 𝜂 1 = 0.01 and 𝜂 2 = 0.9. The whole algorithm is summarized in Algorithm 1.

Advantages of TRRA

• TRRA accept bound constraints

• It uses both gradient and hessian information, making it more robust compared to derivative free algorithms.

• The subproblem 𝑚 𝑘 (𝑠) is defined in lower dimension, thus they are more suitable for large scale compared to other algorithms

Limitations of TRRA

• TRRA may not suitable for all problems, especially problems which are noisy or problems with discontinuous objective function

• The subproblem 𝑚 𝑘 (𝑠) may need to be resolved several times in one iteration before obtaining an acceptable trial step, thus the total cost of computation for one iteration might be expensive.

III.2.1.2 Levenberg-Marquardt (unconstrained)

The Levenberg-Marquardt algorithm (LM) is used to solve non-linear least squares problems.

However, as with many fitting algorithms, the LM finds only a local minimum, which is not necessarily the global minimum. The LM algorithm is a combination of the Gauss-Newton Algorithm (GN) and the Gradient Descent algorithm. The LM algorithm is more robust than the GN algorithm, which means that it will often find a solution even if it starts very far from the minimum. On the other hand, for well-behaved functions and reasonable parameters, LM tends to be a bit slower than the GN. Like other numerical minimization algorithms, the Levenberg-Marquardt algorithm is an iterative procedure used especially for NLS problems. NLS are optimization problems defined as:

min 1 2 ‖𝑟(𝑃)‖ 2 2 ( 24
)
Where 𝑟: ℝ 𝑛 → ℝ 𝑛 is a twice differentiable function.

To solve an NLS problem, the first step is to apply the Newton's method directly. To do this, let us calculate the gradient and the Hessian of the problem. If 𝐶(𝑃) is the cost function defined as: Where 𝐽 𝑟 (𝑃) and 𝐻 𝑟 (𝑃) are respectively the Jacobian and the hessian of 𝑟(𝑃).

𝐶(𝑃) =
This Hessian is in general very expensive to calculate. In addition, the optimum 𝑟 will be small, making the contribution of the term 𝐻 𝑟 (𝑃)𝑟(𝑃) negligible. It is therefore better to ignore it and thus consider that:

𝛻 2 𝐶(𝑃) ⋍ 𝐽 𝑟 (𝑃) ⊤ 𝐽 𝑟 (𝑃)
This approximation is known as the Gauss-Newton hypothesis (defined in chapter 2). It allows on the one hand, to approach 𝛻 2 𝐶(𝑃) by primary derivatives and to ensure, on the other hand, that 𝛻 2 𝐶(𝑃) is always positive semi-definite. This hypothesis leads to the Gauss-Newton algorithm (defined in chapter 2).

However, the convergence of this method is not guaranteed. Indeed, if 𝑃 0 is too far from a local optimum, then 𝐽 𝑟 (𝑃 𝑘 ) ⊤ 𝐽 𝑟 (𝑃 𝑘 ) can be singular or ill-conditioned. Nevertheless, these disadvantages are related to the Newton's method and not to the Gauss-Newton hypothesis itself.

In order to overcome these obstacles, it is therefore possible to solve an approximated quadratic sub-function instead of solving the original cost function. This obstacle could be solved by increasing the diagonal hessian to make it positive and better conditioned.

Advantages of LM

• It uses both gradient and hessian information, which make it more robust compared to derivative free algorithms.

Algorithm 2: Levenberg-Marquardt

Initial conditions: 𝑃 0 ∈ ℝ 𝑛 , 𝜇 0 • LM algorithm combines the advantages of gradient-descent and Gauss-Newton methods.

It inherits the speed advantage of the Gauss-Newton algorithm and the stability of the steepest descent method

• The hessian does not need to be calculated as they are estimated from the gradient of the residuals.

Limitations of LM

• The LM does not handle any type of constraints.

• LM may not be suitable for all problems, especially problems which are noisy or problems with discontinuous objective function

III.2.1.3 Nelder-Mead Algorithm

The Nelder-Mead algorithm (NM) or simplex search algorithm, originally published in 1965 [5],

is one of the best known algorithms for multidimensional unconstrained optimization without derivatives. The NM is quite simple to understand and implement in practice, because it does not require gradient computation (i.e., derivative free), which makes it suitable for problems with nonsmooth functions or blackbox optimization problems. For these reasons, it is very popular in many fields of science and technology, especially in chemistry, physics, and medicine [6]- [8]. It is widely used to solve parameter estimation and statistical problems where the function values are uncertain or subject to noise. It can also be used for problems with discontinuous functions, which occurs frequently in statistics and experimental mathematics [9].

Assume 𝐶(𝑃): ℝ 𝑛 → ℝ be the cost function to be minimized. This method starts with a simplex (i.e., polytope) of 𝑛 + 1 vertices, which is a point of ℝ 𝑛 , denoted 𝑃 0 , 𝑃 1 … , 𝑃 𝑛 , such that 𝐶(𝑃 0 ) ≤ 𝐶(𝑃 1 ) ≤ ⋯ ≤ 𝐶(𝑃 𝑛 ) , see Figure 17.

We are trying to minimize 𝐶, 𝑃 0 is the best point and 𝑃 𝑛 is the bad one. At each iteration, the bad point 𝑃 𝑛 is eliminated, and another point is accepted into the simplex. This process continues until convergence is achieved. An example of the process is illustrated in Figure 17, 18, 19, 20 and 21.

Suppose that our simplex the one represented in Figure 17, such that: 𝐶(𝑃 0 ) ≤ 𝐶(𝑃 1 ) ≤ 𝐶(𝑃 2 )

Chapter The algorithm is detailed in the Algorithm 3.

Algorithm 3: Nelder-Mead

Initial conditions: 𝑛 + 1 points 𝑃 0 , … , 𝑃 𝑛

Advantages of Nelder-Mead

• Does not require any derivative information (i.e., derivative free problems).

• NM is quite simple to understand and very easy to implement in practice

• Optimizes both continuous and discrete functions and also suits multi-objective problems.

• The Nelder-Mead method frequently gives significant improvements in the first few iterations and quickly produces quite satisfactory results.

• The method typically requires only one or two function evaluations per iteration (i.e., This is very important in applications where each function evaluation is very expensive or timeconsuming)

Limitations of Nelder-Mead

• The method can take an enormous number of iterations with negligible improvement in function value.

• Extremely depending on the choice of the optimization starting point

III.2.2 Stochastic Algorithms

Stochastic optimization algorithms are algorithms that make use of randomness in the search procedure for objective functions for which derivatives cannot be calculated. Many stochastic algorithms are inspired by biological or natural processes and may be referred to as "metaheuristics" as a higher-order procedure providing the conditions for a specific search of the objective function. They are also referred to as "black box" optimization algorithms.

Stochastic optimization algorithms provide an alternative approach that allows less optimal local solutions to be found within the search procedure, which may increase the probability of the procedure locating the global optima of the objective function. Algorithms of this type are proposed for the most challenging optimization problems that may have noisy function evaluations or many global optima (multimodal). Two kinds of stochastic optimization algorithms will be addressed in this chapter, which are, Genetic Algorithm and Particle Swarm Algorithm.

III.2.2.1 Genetic Algorithm

The Genetic Algorithm (GA) [10] is a global search heuristic technique based on the principles of genetics and natural selection. It is frequently used to find optimal solutions to difficult problems that may take a lifetime to solve. It employs techniques that are inspired by evolution, such as mutation, crossover, and selection. The algorithm searches for the best solution in the search space that optimizes the cost function. GA does not require gradient computation of the cost function (i.e., derivative free), which makes it suitable for problems with non-smooth functions or blackbox optimization problems. The algorithm starts with randomly generated parameters, then the parameters are tested using a fitness value which is obtained from the cost function, and finally the best ones are selected probabilistically from the current iteration. Iteration terminates when the desired fitness value is obtained or the maximum number of iterations has been reached.

Briefly, the iterative procedure of GA is as following: The key point of this algorithm is the step where the new parameters are reproduced by crossover and mutation. The iterative procedure of the genetic algorithm's evolution is illustrated in Figure 22.

Advantages of Genetic Algorithm

• Does not require any derivative information (i.e., derivative free problems).

• More accurate than many traditional methods.

• Optimizes both continuous and discrete functions and also multi-objective problems.

• Suitable for problems with a large search space or many parameters.

Limitations of Genetic Algorithm

• GAs are not suitable for all problems, especially problems which are simple and for which derivative information is available.

• Cost function is calculated repeatedly which might be computationally expensive for some problems.

• If not implemented properly, the GA may not converge to the optimal solution.

Chapter

III.2.2.2 Particle Swarm Algorithm

Particle Swarm algorithm (PS) [11] has become one of the most promising optimization techniques for solving global optimization problems. In the past few years, PS has been successfully applied in many research and application areas [12], [13]. PS is a metaheuristic as it makes few or no assumptions about the problem being optimized and performs very well on a wide range of problems, including high-dimensional problems and some that are very difficult to solve using classic numerical optimization techniques. Also, PS does not require gradient computation of the cost function (i.e., derivative free), which makes it suitable for problems with non-smooth functions or blackbox problems. However, metaheuristics such as PS do not guarantee that a global optimal solution can be found for some class of problems.

The PS system solves problems by having a population (swarms) of candidate solutions (particles).

These particles move through the search domain with a specified velocity in the search for an optimal solution. Each particle keeps track of its coordinates in the search space, which are associated with the best solution it has achieved so far. This value is called "personal best." Another value that is tracked by the PS is the best value obtained so far by any particle in the neighbors of Chapter III 90 the particle. This location is called local best. When a particle takes all the population as its topological neighbors, the best value is called global best.

Let 𝐶(𝑃): ℝ 𝑛 → ℝ be the cost function to be minimized with 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏 . The function takes a candidate solution as an argument in the form of a vector of real numbers and produces a real number as output which indicates the objective function value of the given candidate solution. The goal is to find a solution 𝑃 ⋆ for which 𝐶(𝑃 ⋆ ) ≤ 𝐶(𝑃) for all 𝑃 in the search-space.

Let 𝑛 be the number of particles in the swarm, each having a position 𝑥 𝑖 ∈ ℝ 𝑛 in the search-space and a velocity v i ∈ ℝ 𝑛 . At iteration 𝑘, the position of particle 𝑖 = 1, … , 𝑛 is denoted by 𝑥 𝑖 (𝑘) and the velocity by v 𝑖 (𝑘). The position of particle 𝑖 at iteration 𝑡 is simply the previous position (at iteration 𝑘 -1) plus the velocity at the previous iteration, such that:

𝑥 𝑖 (𝑘) = 𝑥 𝑖 (𝑘 -1) + v 𝑖 (𝑘 -1)
The fitness of a position is determined by the value of the cost function at that position. The best position visited by a particle, 𝒑 𝒊 (𝒌), is simply the position visited by particle 𝑖 in iterations 1, … , 𝑘 with the minimal function value, i.e., it is the position corresponding to the optimal value of 𝑝 𝑖 (𝑘) = min{𝐶(𝑥 𝑖 (0)), … , 𝐶(𝑥 𝑖 (𝑘))}, where 𝐶 is the cost function. Similarly, the best position visited by any particle in the swarm, known as the global best, 𝒈(𝒌), is the position corresponding to the optimal value of:

{𝐶(𝑥 1 (0)), … , 𝐶(𝑥 1 (𝑘)), 𝐶(𝑥 𝑛 (0)), … , 𝐶(𝑥 𝑛 (𝑘))}.

Thus, the velocity of particle 𝑖 at iteration 𝑡 is:

𝑣 𝑖 (𝑘) = 𝑐 0 𝑣 𝑖 (𝑘 -1) + 𝑐 1 𝑟 1 (𝑥 𝑖 (𝑘) -𝑝 𝑖 (𝑘)) + 𝑐 2 𝑟 2 (𝑥 𝑖 (𝑘) -𝑔(𝑘))

Where 𝑟 1 and 𝑟 2 are chosen randomly by a uniform distribution on the interval [0,1]. The constants 𝑐 0 , 𝑐 1 and 𝑐 2 weight each of the three components of the particle's velocity, which are defined to control the behavior and efficiency of the PS method.

PS iterative procedure is detailed in Algorithm 5. The procedure of updating particle position and velocity is illustrated in Figure 23.

Chapter

Advantages of Particle Swarm

• Does not require any derivative information (i.e., derivative free).

• A very efficient global search algorithm.

• PS is less sensitive to the nature of the objective function compared to the predictable mathematical approaches and other heuristic methods

• PS has limited number of parameters including only inertia weight factor and two acceleration coefficients in comparison with other competing heuristic optimization methods.

• PS techniques can generate high-quality solutions within shorter calculation time and stable convergence characteristics than other stochastic methods [14].

Limitations of Particle Swarm

• It requires a longer computation time compared to direct search algorithms.

III.3 Comparison between optimization algorithms

It is worth mentioning that the examples shown in the following sections are a simplified version of our model (i.e., only 8 unknown parameters are considered and only a small period of polarization time is considered, see Table 3). When the optimization technique is used on a large amount of data, the computation time may become prohibitive. That is to say, one of the important aspects that should be investigated to determine algorithm efficiency is time consumption.

Subsequently, several initial points have been tested to validate the optimization outcomes, but only one of them has been chosen and presented in this section. The example was chosen in which it might summarize all the possibilities we could face during the optimization process. The upper and lower bounds of the parameters (Table 3) are considered based on the experimental approaches made by various researchers [15]- [17]. In this section, only the charge density data will be studied since it provides more information than the current density data.

The methodology is as follows, assume that we want to find the value of a set of parameters vector of real variables 𝑃 ⋆ such that it minimizes the cost function 𝐶(𝑃), defined as:

𝐶(𝑃) = ‖𝑛 𝐸𝑥𝑝 -𝑛 𝑆𝑖𝑚 (𝑃)‖ (25)
Where 𝑛 𝑆𝑖𝑚 (P) are the simulated charge density by the BCT model in terms of unknown parameters P and 𝑛 𝐸𝑥𝑝 represents the experimental charge density.

The optimization problem to be solved is thus defined as: min 𝑃∈ℝ 𝑛 𝐶(𝑃) 𝑠. 𝑡. 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏 [START_REF] Zerbi | Transport of matter upon annealing in crystalline polymethylene systems: a calorimetric and spectroscopic study[END_REF] Chapter III 94 Five different optimization algorithms are used to solve Problem [START_REF] Zerbi | Transport of matter upon annealing in crystalline polymethylene systems: a calorimetric and spectroscopic study[END_REF]. All the algorithms were implemented and coded in MATLAB. TR and LM were implemented using the "lsqnonlin" function. GA was implemented using the "ga" function. NM algorithm was applied using the "fminsearchbnd" function. Finally, the PS algorithm was applied using the "particleswarm" function.

III.3.1 Optimization using simulated data

The major purpose of this section is to examine five alternative optimization algorithms to validate and evaluate their convergence performance by applying them on the BCT model, then compare their outputs to determine the most appropriate one for solving our problem.

In this part, the experimental data are simulated using a bipolar charge transport model using a set of parameters known as 𝑃 𝐸𝑥𝑎𝑐𝑡 (displayed in Table 2), which are obtained from the literature [18]. Thus, the cost function will be defined as:

𝐶(𝑃) = ‖𝑛 𝑆𝑖𝑚 (𝑃 𝑒𝑥𝑎𝑐𝑡 ) -𝑛 𝑆𝑖𝑚 (𝑃)‖ (27)
Where 𝑛 𝑆𝑖𝑚 (P) are the simulated charge density by the BCT model in terms of unknown parameters P and 𝑛 𝑆𝑖𝑚 (𝑃 𝑒𝑥𝑎𝑐𝑡 ) represents the charge density data that are obtained by simulating the BCT model using simulated by 𝑃 𝑒𝑥𝑎𝑐𝑡 .

The model in this part intends to model the behavior of 200 µ𝑚 LDPE under a constant DC stress of -30kV/mm for 500 s of polarization time at 25 °𝐶 using a set of parameters 𝑃 𝐸𝑥𝑎𝑐𝑡 , that are displayed in Table 4. Table 4 represents the optimal set of parameters produced by each one of the used optimization algorithms. 𝑃 𝑇𝑅𝑅𝐴 , 𝑃 𝐿𝑀 , 𝑃 𝑁𝑀 , 𝑃 𝐺𝐴 , and 𝑃 𝑃𝑆 , are the optimal parameters produced respectively by, TRRA, LM, NM, GA, and PS. Table 5 represents the time of computation, the number of iterations, and the number of evaluated functions required by each algorithm to find the optimal solution. Also, it represents the relative error which evaluates the cost function using the obtained optimal solutions. Based on Table 4 and Table 5, GA and PS took a very long time to converge toward the optimal solution, compares to other algorithms. Besides, PS produced a very high relative error (i.e., the convergence of 𝑃𝑆 was found to be bad in comparison to other algorithms, where 𝐶(𝑃 𝑃𝑆 ) = 43.68 is very high compared to the other algorithms). According to Table 4, comparing 𝑃 𝑃𝑆 with 𝑃 𝐸𝑥𝑎𝑐𝑡 proofs the bad convergence of PS algorithm, where most of the parameters did not converge to the exact value.

Chapter III 95

Besides, based on Table 4, it is obvious that TRRA is the most accurate algorithm compared to the others. This could be noticed by comparing the relative error produced by each optimal solution (𝐶(𝑃 𝑇𝑅𝑅𝐴 ) < 𝐶(𝑃 𝐿𝑀 ) < 𝐶(𝑃 𝐺𝐴 ) < 𝐶(𝑃 𝑁𝑀 ) < 𝐶(𝑃 𝑃𝑆 )). This also could be noticed from Table 4, where all the parameters achieved their exact value when TRRA was used, while for other algorithms, most of the parameters did not converge to the exact value.

Furthermore, based on Table 5, its clear that LM has the lowest computation time and number of iterations, compared to other algorithms. Also, according to the relative error produced by 𝑃 𝐿𝑀 Chapter III

96

(𝐶(𝑃 𝐿𝑀 )), it is found that 𝑃 𝐿𝑀 is near to 𝑃 𝐸𝑥𝑎𝑐𝑡 , where most of the parameters converged to their exact values.

Finally, the parameters produced by NM (𝑃 𝑁𝑀 ) are not compatible with the exact parameters (𝑃 𝐸𝑥𝑎𝑐𝑡 ). Although, when the NM method is performed using a different starting point (𝑃 0 ), the outcomes are identical to 𝑃 𝐸𝑥𝑎𝑐𝑡 . That is to say, the NM method was found to be extremely dependent on the starting point of the optimization processes.

To summarize what has been stated so far, TRRA and LM was found to be the most appropriate algorithms for solving our problem. TRRA had the most accurate optimal solution, whereas LM converged faster than all the other algorithms. GA produced an acceptable optimal solution but it required a lot of computation time to be achieved. PS produced a very weak optimal solution with a lot of time to be achieved. Finally, NM converged quickly but the produced optimal parameters did not match 𝑃 𝐸𝑥𝑎𝑐𝑡 , since NM was found to be very sensitive with respect to the choice of initial starting points. Table 6 shows the outcomes of the NM algorithm using two different starting points. The results displayed in Table 6 show that if the initial starting point 𝑃 0 was chosen near to the optimal solution (Example 1), NM can find the exact value of most of the parameters.

However, if 𝑃 0 was chosen randomly (Example 2), NM may not achieve the exact values of the parameters we are seeking to find.

III.3.2 Optimization using experimental data (PEA)

Another approach is used in this subsection to compare and validate the optimization algorithms so that the best one will be chosen. This approach aims to optimize the same cost functions as in the previous subsection, but with PEA experimental data rather than model data.

Thus, the cost function will be defined as: The outcomes of the optimization algorithms are displayed in Table 7 and Table 8. Table 7 represents the optimal set of parameters produced by each one of the used optimization algorithms. 𝑃 𝑇𝑅𝑅𝐴 , 𝑃 𝐿𝑀 , 𝑃 𝑁𝑀 , 𝑃 𝐺𝐴 , and 𝑃 𝑃𝑆 , are the optimal parameters produced respectively by, TRRA, LM, NM, GA, and PS.

𝐶(𝑃) = ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑆𝑖𝑚 (𝑃)‖ (28)
Table 8 represents the time of computation, the number of iterations, and the number of evaluated functions required by each algorithm to find the optimal solution. Also, it represents the relative error which evaluates the cost function using the obtained optimal solutions.

Chapter III 98 Based on Table 7, it is obvious that TRRA is the most accurate algorithm compared to the others.

This could be noticed by comparing the relative error produced by each optimal solution (i.e., 𝐶(𝑃 𝑇𝑅𝑅𝐴 ) < 𝐶(𝑃 𝐿𝑀 ) < 𝐶(𝑃 𝐺𝐴 ) < 𝐶(𝑃 𝑁𝑀 ) < 𝐶(𝑃 𝑃𝑆 )). Besides, TRRA is found to be the most efficient algorithm compared to the others, where the optimal solution was found with minimal computation time compared to the other algorithms. The findings in Table 7 reveal that LM is not suitable for optimizing the BCT model, because the value of holes trapping coefficient (𝐵 ℎ ) was found with a negative sign, which is wrong because the trapping coefficient is known physically as a positive parameter. This is due to the fact that the LM algorithm does not handle bound constraints. This makes LM a bad choice for solving our problems.

Besides, it could be seen from Table 8 that the LM consumed more time than TRRA, even that the number of iteration (=26 iterations) needed by LM is less than that of needed by TRRA (=57 iterations). This is due to the absence of the bounds for LM algorithm, and then the parameters converged a bad optimal solution that slows the BCT model performance at each iteration.

Furthermore,

Table 8 shows that GA and PS produced an accurate optimal solution based on the value of the relative error (i.e., 𝐶(𝑃 𝐺𝐴 ) = 0.133 and 𝐶(𝑃 𝑃𝑆 ) = 0.101)) but the computation time was prohibitive for both algorithms.

NM algorithm produced an accurate optimal solution (𝐶(𝑃 𝑁𝑀 ) = 0.111) with an acceptable computational time (=3945 s), however, the NM failed to converge into the exact value of the parameters in the previous subsection. That is to say, NM algorithm extremely dependent of the initial starting point. Whereas, TRRA was always able to find an acceptable set of parameters with any initial set of parameters. This shows that TRRA is more robust than NM algorithm while optimizing the BCT model. simulated by using P TRRA , LM: simulated by using P LM , NM: simulated using P NM , GA: simulated by using P GA , and PS: simulated by using P PS . According to Figure 25, it is obvious that the simulated results obtained by 𝑃 𝑇𝑅𝑅𝐴 and 𝑃 𝑃𝑆 provide the best correlation between experimental and simulated cartography's. These results are comparable with results presented in Table 8, where the minimal relative error was found for 𝑃 𝑇𝑅𝑅𝐴 (𝐶(𝑃 𝑇𝑅𝑅𝐴 ) = 0.083) and 𝑃 𝑃𝑆 (𝐶(𝑃 𝑃𝑆 ) = 0.101).

III.4 Conclusion

The parameters of the charge transport model were optimized using five different optimization algorithms. In order to select the most reliable algorithm, a trivial example was employed to examinate all the considered algorithms. To find the best algorithm, two distinct approaches were employed. The first approach optimizes parameters using modeled data, whereas the second approach optimizes experimental charge density data obtained by using the PEA method. TRRA was found to be the most suitable algorithm for this type of problem in both approaches (Table 9).

When compared to the other methods, TRRA had the best convergence rate and the shortest calculation time. NM might also be a good option, but it is totally dependent on the initial point selection, making TRRA the preferable option. The accuracy of GA is acceptable; however, the calculation time is quite long. LM was found to be efficient but it may converge to an unsatisfactory solution due to the fact that the LM method does not handle bound constraints.

✓ ✓ ✓ LM x ✓ x NM x ✓ x GA ✓ x x PS x x x

IV.1 Introduction

In the previous chapter, several optimization algorithms were applied to a simple example of our problem to compare them and choose the best one. Trust Region Reflective Algorithm (TRRA) was found to be the most suitable algorithm to solve our problem and achieve our target. Thus, the main target of this chapter is to apply TRRA considering all the unknown parameters we are seeking to find.

Indeed, the proposed model which is described in chapter 1, requires some experimental conditions such as temperature and applied voltage, alongside a set of unknown parameters, including injection barrier, mobility, trapping and detrapping rates, and recombination coefficients. Most of these parameters cannot be predicted, detected, or estimated using independent experiments. For this reason, this chapter proposes an original method based on a dual approach (Simulation/Experiment) that helps to critically analyze the BCT model (i.e., highlighting the weaknesses and strengths) and to have a good correlation between simulated and experimental behavior. This could be done by minimizing the difference between experimental and simulated data for both current and charge experiments using TRRA with several electrical fields.

This study considers the two most accessible experiments in our laboratory, which are the net charge density that is obtained using the PEA method, and the current density obtained from external charging current measurements. Many other experimental sources could be considered in future to enhance our approach, such as electroluminescence. On the other side, the BCT model established for LDPE under DC stress is used to estimate the simulated charge and current densities.

This study shows the importance and influence of incorporating current and charge density measurements together with varying electric fields on the optimization algorithm behavior. For simplicity reasons, the majority of the proposed models in the literature considered that the trapping and detrapping rates remain constant as the electric field increases [1]- [3]. However, experiments have widely established that these rates are proportional to the electric field [4], [5].

For this reason, TRRA is applied in this paper using several experimental measurements with varying applied fields to understand the impact of the applied electric field on the trapping and detrapping rates.

Chapter IV

106
This chapter is composed of five sections. After this brief introduction, the second section presents the sample conditioning and the experimental set ups utilized for this study, with an example illustrating the outputs of each experiment. In the third section, we will present the strategy used to optimize the BCT model based on specific experimental protocol. In the fourth section, TRRA is applied using several cost functions to find the optimal set of parameters that fit all of the provided experiments. Finally, the results are discussed and analyzed in depth.

IV.2 Sample preparation and characterization tools IV.2.1 Material preparation

The LDPE samples are manufactured in the laboratory using a heating press with BOREALIS polyethylene resin pellets. Since our study aims at performing all experiments using the same LDPE film, it was thus manufactured with specific characteristics to be suitable for all used experiments (i.e., the LDPE film thickness should be in the range [100 μm -200 μm]). The thickness of the samples depends on the mass of pellets, the temperature, and the pressure during the manufacturing protocol.

The LDPE pellets are sandwiched between two polyimide films that are placed on the aluminum disc (2-3 inches) of the heating press (Figure 26). These films prevent the LDPE sample from a direct contact with the press aluminum disc, which are not smooth. They also allow for easier removal of the LDPE sample after processing. A 200 μm Polytetrafluoroethylene (PTFE) spacer with an internal diameter of 80 mm is placed between the polyimide films and serves as a mold for the LDPE. Figure 27 shows the pressure and temperature cycle used to prepare the LDPE samples. No pressure is applied on the pellets until reaching a temperature of 155 °C (10 minutes), in order to melt them. The pressure is then progressively increased to 2.5 tons, while remaining at 155 °C for 10 minutes. Finally, a temperature set point of 30 °C is set. It takes around 40 minutes for the cooling process to complete. To avoid the formation of bubbles, the pressure remained steady at 2.5 tons until the temperature reached 30 °C. For preparing the LDPE sample, we used 0.8 g of LDPE pellets to obtain disks of 150 μm thickness, and 60 mm diameter (Figure 28). For current measurements, the LDPE should be metallized by coating the sample by a surface electrode on both sides. The process of metallization is presented in the following section.

Chapter

IV.2.2 Sputtering metallization

Sputtering metallization technique was first observed in 1852, during the operation of discharge tubes, and since then, this method has been widely used for thin film metal deposition [6]. The mechanism of sputtering is based on moment transfer [7]. A DC voltage is applied between the metal target (cathode) and the substrate (anode) after evacuation of air and introduction of argon at a pressure of about 6 × 10 -2 mbar. The cathode-anode gap can vary according to the models from 4 to 10 cm. The anode holder is cooled by water circulation. Under the effect of the electric field, the ions are attracted by the cathode (target) and the electrons move towards the anode and maintain the plasma. Figure 29 shows the sputtering device used for the gold deposits.

The metallization is essential for current measurements to ensure a good contact between the electrodes and LDPE film (i.e., current measurements are very sensitive and require a good contact between the sample and electrodes). Furthermore, for space charge measurements, the LDPE was metallized to provide symmetrical electrodes, the same as for current measurements (i.e., the upper electrode of a PEA cell is a semiconductor and the lower electrode is aluminum).

In this study, the gold deposit was performed by sputtering on both sides of the LDPE samples.

The LDPE sample is sandwiched between two masks that are used to choose the shape of the gold deposit, which is in our case a 10 mm diameter circle (see Figure 30). With this method, the thickness of the deposit layer is not homogeneous in the plane of the sample and follows the shape of a Gaussian: thin at the edges (30 nm) and thicker in the center (50 nm) of the sample.

IV.2.3 Current measurement

The cell used for current measurements consists of a flat sample of insulating material of thickness D that is sandwiched between two electrodes with 1 cm diameter that have been optically polished (see Figure 31). A positive or negative DC voltage 𝑉 0 is applied to the upper electrode at zero time, while the other electrode is grounded. Furthermore, Keithley 6512 is used to measure the current 𝐼(𝑡) during the application of the electrical stress.

𝐽 𝑐 (𝑡) = 𝐼(𝑡) 𝑆 ( 29
)
Where S is the area of the surface electrode.

The curve displayed in Figure 32 According to Adamec and Calderwood [8], a weak polarization process is dominant in LDPE at short periods and it is a function of the applied field. The polarization of the dipoles inside the bulk might explain the initial slope at short time (𝑡 < 10 seconds). Moreover, at short time, the current measurements are unreliable due to the used electrometer. For these reasons, the first 10 seconds won't be considered.

IV.2.4 Pulsed Electro-Acoustic method: Charge density measurement IV.2.4.1 Principle of PEA method

The presence of space charges in materials is a phenomenon which has been studied for many years [9], [10]. However, understanding of the generation, transport, and trapping modes of these charges has not yet been achieved. Amongst the numerous techniques developed to detect the presence of space charges in materials [11], this chapter focuses on the PEA method. It is an acoustic technique for the detection of the net density of charge (i.e., the algebraic sum of all electrical charges). It was developed in Japan by Maeno and his colleagues in 1988 [12]. The sample is placed between the two electrodes of the PEA cell, the lower electrode made of aluminum and the upper electrode is made of semi-conductor. The PEA cell could be placed in a thermal chamber in order to control the measurement temperature. The principle of this method consists in applying electrical pulses with a pulse generator to a sample placed under DC or AC voltage. Under the effect of the pulses, the space charges in the dielectric move around their point of origin. Figure 33 shows the principle of PEA for a negative charge implanted in the sample.

These negative charges induce a positive image charges on each electrode. When a pulse electric field up(t) is applied across the sample, an acoustic wave is created by displacement of the charges around their original position. These signals contain information about the position of the charges in the volume of the dielectric as well as their quantity. This wave 𝑃 ∆ (t) propagates at a speed of Chapter IV 111 sound in the material and is collected by a piezoelectric sensor. This piezoelectric sensor translates the wave 𝑃 ∆ (t) into an electrical signal vs(t). Then, a signal processing is performed to estimate the space charge profiles as a function of the thickness and time.

IV.2.4.2 Experimental device

The schematic diagram of the PEA device is shown in Figure 34. It is made up of the following elements:

1) A pulse generator that provides an electrical signal with a frequency of 1kHz and amplitude 250V in order not to influence the charge distribution in the material.

2) A high voltage DC source used to polarize the sample.

3) A digital oscilloscope and a computer for data acquisition and digital processing, respectively.

The upper electrode of the PEA cell is made of semi-conductor, with a diameter 10 mm and the lower electrode is made of aluminum. Below the lower electrode, there is a 9 μm PVDF (polyvinylidene fluoride) piezoelectric sensor. The PEA cell we use has a spatial resolution of about 10 μm and a sensitivity of 0.1 C/m 3 . All the space charge measurements presented in this manuscript have been performed using this method. An example of a charge density measurement using PEA is represented in Figure 35.

Chapter that is [-2 , 2] 𝐶/𝑚 3 . In this example, the cathode is the upper electrode and the anode is the lower electrode, such that electrons are injected from the cathode and holes are injected from the anode.

According to Figure 35, the positive and negative charges are detected immediately after the voltage is applied. The dominance negative charge could be seen with a powerful and quick

Chapter IV 113 injection over a short period of time (Zone 1). After 600 seconds (Zone 2), the positive charges have the advantage inside the bulk of insulation, compared to negative charges.

IV.3 Experimental results

IV.3.1 Experimental protocol

In this study, two different experiments have been carried out using the same LDPE sample.

Experiment 1 aims to measure the charge density in LDPE at three levels of fields of -20, -40, and -60 𝑘𝑉/𝑚𝑚. The same protocol is used for Experiment 2 for measuring the current density.

Both experiments consist of polarization of 3 hours and a short-circuiting period of 24 hours each, as shown in Figure 36. The sample has been discharged 24 hours after each measurement to eliminate the charges in the bulk as much as possible. Space charge measurements were performed using the PEA technique. The current density was obtained by the known external current measurements. The experimental protocol is illustrated in Figure 36. Since the study is based on a dual approach (experiment/simulation), it was better to perform all the measurements using the same sample and exactly on the same location. This was done to ensure that all experiments were done using the same thickness and same material properties. Since all the measurements were performed on the same sample, the sample was discharged for 24 hours after each experiment to remove as many charges as possible from the bulk. However, it is well known that it is difficult to eliminate all the charges from the material. For this reason, a finite charge density can be set as an initial condition in the simulated model (𝛼) with 𝐶. 𝑚 -3 as a unit.

This level of charge is assumed to simulate a residual charge which is likely to exist in any Chapter IV 114 dielectric but is not seen in space charge measurements due to the limited sensitivity of the detection techniques (order of 0.1 𝐶. 𝑚 -3 ) or to the fact that densities of positive and negative charges are locally equal, giving a zero net charge density. Regardless of the applied voltage, the dominance negative charge could be seen with a powerful and quick injection over a short period of time. After 500 seconds, the dominance of positive charges could be observed inside the bulk of insulation. Moreover, the mobility of electron seems to be higher than that of holes at any electric field.

IV.3.2 Experimental measurements

Figure 38 represents the net density of charge for -60 𝑘𝑉/𝑚𝑚 applied electric field at 2 periods of time, 100 s (blue) and 1000 s (red). Based on Figure 38, it is obvious that the amount of the negative charges is more than the amount of positive in the bulk of LDPE at short time (𝑡 < 500 𝑠). After 500 seconds, the amount of positive charges increased inside the bulk of insulation, whereas the amount of negative charges decreased. As mentioned previously, the first 10 seconds of experimental current measurements will be ignored during the optimization process. represents the initial net charge density before starting current measurements for fields -20, -40 𝑎𝑛𝑑 -60 𝑘𝑉/𝑚𝑚, respectively. The units, symbols, and bounds of all the unknown parameters are displayed in Table 10. The lower and upper bounds are chosen based on previous experimental measurements made by various researchers, such that:

• For electron injection, the theoretical injection barrier height corresponding to a goldpolyethylene interface is 5 eV. When this value was applied to the Schottky law, the 3 eV] based on the outcomes of charge transport models in polyethylene [13], [14]. Therefore, this range will be considered in our study.

• The PEA measurements have been used by several researchers [15], [16] to approximate the charge mobility. Depending on the type of carrier and the temperature, the final estimation produced mobilities in the range of 10 -12 -10 -14 𝑚 2 𝑉 -1 𝑠 -1 . Thus, the ideal range of shallow trapping depth is [0.3 -0.72 𝑒𝑉] and inter-trap distance in the range

[0 -20 𝑛𝑚].
• Some of the parameters, such as recombination coefficients, could not be defined easily by experiments. For this purpose, the criticality of the parameters was examined, i.e., by studying the effects of modifying the given parameter on the model outputs [17]. It has been shown that the recombination coefficients do not play a major role in simulating results at low temperatures or short periods of polarization. By considering one of these cases, recombination coefficients do not change the current density to a significant degree, and they weakly affect the charge density in the dielectric, in the limit of the investigated range of recombination rates from 10 -5 to 1 m 3 /Cs. Thus, the ideal range of recombination pre-factors is

[10 -3 -1].
Furthermore, in the literature, almost all of the present models consider the trapping and detrapping rates constant with respect to the increase of the electric field [1]- [3]. Nevertheless, based on experimental studies, it is suspected that these rates are related to the electric field, temperature, or other experimental conditions [4], [5]. It is worth mentioning that the detrapping phenomenon has previously been described using an equation that estimates the rate of detrapping as a function of temperature. However, no equations exist that could describe trapping and detrapping phenomena as a function of the electric field.

It is therefore essential to better understand trapping and detrapping processes in the material. For this reason, twelve additional parameters were introduced into the optimization process, where each one represents the trapping and detrapping coefficients at a different electric field. The additional parameters are defined in Table 11.

Where 𝐵 𝑒 𝑖 , 𝐵 ℎ 𝑖 are the trapping rates of electrons and holes, respectively, for 𝑖 = -20, -40 and -60 𝑘𝑉/𝑚𝑚 represents the applied fields used for our measurements. 𝐷 𝑒 𝑖 , 𝐷 ℎ 𝑖 are the detrapping rates of electrons and holes, respectively, for 𝑖 = -20, -40 and -60 𝑘𝑉/𝑚𝑚 represents the applied fields used for our measurements. The intensity of the deep trapping process is described by the detrapping coefficient 𝐷 𝑒,ℎ , which reflect the detrapping rate per unit of time, are defined as:

𝐷 𝑡𝑟𝑒 = 𝐷 𝑒 ⋅ 𝑛 𝑒𝑡 (30)
𝐷 𝑡𝑟ℎ = 𝐷 ℎ ⋅ 𝑛 ℎ𝑡 (31)
Previous studies proposed an equation to estimate the detrapping rates 𝐷 𝑒,ℎ with respect to temperature, such that:

𝐷 𝑒,ℎ = 𝑣 • 𝑒xp ( -𝑤 𝑡𝑟 𝑒,ℎ 𝑘 𝐵 𝑇 ) (32)
Where 𝑣 is the frequency of the phonons, 𝑤 𝑡𝑟 𝑒,ℎ are the barrier heights to escape from deep traps, 𝑘 B = 1.381 × 10 -23 𝐽/𝐾 stands for Boltzmann constant, and T is the temperature.

In this study, instead of using Equation 5, the detrapping rates 𝐷 𝑒,ℎ will be supposed as an unknown parameter that varies depending on the electric field (Table 11), because the temperature is constant.

In total, 27 unknown model parameters are included in the optimization process. The trapping and detrapping phenomena are described by using 12 distinct parameters (as shown in Table 11), to evaluate their values at different electric fields. 6 additional parameters are introduced to the optimization process to evaluate the initial charge density at the beginning of each experiment.

Besides, 2 parameters are used to represent the injection barrier heights of electrons and holes. 2

Chapter IV 119 parameters represent the shallow trap depth of electrons and holes. 2 represent the inter-trap distance between two shallow traps. Finally, 3 parameters represent the recombination process between different types of carriers.

In our case, TRRA is used to estimate the unknown variables to achieve the minimum squared difference between the experimental and simulated data (Figure 41). It was implemented and coded in MATLAB using "lsqnonlin" function. The first and most important step to start an optimization method is to select the starting point of the optimization algorithm. A bad choice of the starting point, especially in our case (i.e., large number of unknown parameters), may produce a prohibitive time of computation, as well as the convergence rate that may be also affected (Figure 42). In order to avoid this, a particular technique has been used to produce an initial approximation of 𝑃 0 that has the potential to converge faster to the optimal solution. After studying the sensitivity of each parameter with respect to time using Sobol's method [19], results show that exact values of 𝑤 𝑒 and 𝑤 ℎ can be well approximated by optimizing the data of the first [10 -300] seconds of polarization (i.e., the injection parameters are the most effective parameters at short time). In other words, optimizing the difference between experimental and simulated results, using the data of the first 300 seconds, will produce a good approximation for 𝑤 𝑒 and 𝑤 ℎ , and thus, a good approximation for the starting point 𝑃 0 . For this reason, the optimization process will be divided into two steps. First, all parameters will be optimized, considering only the first 300 seconds. Then, the optimal set of parameters produced by the first step will be used as the starting point to optimizing the data considering the total polarization time.

IV.4.2 Strategy to choose a cost function

The methodology is as follows: assume that we need to find the value of a set of parameters vector of real variables 𝑃 ⋆ such that it minimizes the cost function 𝐶(𝑃), defined as:

𝐶(𝑃) = ‖𝑓 𝐸𝑥𝑝 -𝑓 𝑆𝑖𝑚 (𝑃)‖ (33)
Where 𝑓 𝑆𝑖𝑚 (𝑃) are the simulated data by the BCT model and 𝑓 𝐸𝑥𝑝 represents the experimental data. 𝑓 𝐸𝑥𝑝 can be substituted by several experimental measurements, such as, charge density, current density, or electroluminescence, etc.. In this study, we will only focus on using the charge and current density measurements.

Thus, the optimization problem to be solved is defined as:

min 𝑃∈ℝ 𝑛 𝐶(𝑃) 𝑠. 𝑡. 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏 (34)
The major objective of this research is to identify a unique set of parameters that can fit both current and charge experiments, regardless of the used experimental conditions (EC).

It has recently been proved [20] that integrating multiple experiments with different experimental conditions and applying them in the optimization process can help in enhancing the convergence accuracy and efficiency of the optimization algorithm toward the globally optimal set of parameters (Figure 43). Therefore, in this section the TRRA will be used to optimize a multi-objective cost function. This cost function can be formulated in 3 different ways, such that: Case 1: The cost function uses only the data provided by the current density measurements. min 𝑃∈ℝ 𝑛 𝐶 1 (𝑃) 𝑠. 𝑡. 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏 [START_REF] Roy | Modelling space charge in a cable geometry[END_REF] Chapter IV 122 where

𝐶 1 (𝑃) = ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -20𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -20𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -60𝑘𝑉⋅𝑚𝑚 -1
Case 2: The cost function uses only the data provided by the charge density measurements.

min 𝑃∈ℝ 𝑛 𝐶 2 (𝑃) 𝑠. 𝑡. 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏 (36)
Where

𝐶 2 (𝑃) = ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -20𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -20𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -60𝑘𝑉⋅𝑚𝑚 -1
Case 3: The cost function uses the data for both charge and current density measurements (i.e., the combination of case 1 and case 2)

min 𝑃∈ℝ 𝑛 𝐶 3 (𝑃) 𝑠. 𝑡. 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏 (37)
Where:

• 𝐶 3 (𝑃) = 𝐶 1 (𝑃) + 𝐶 2 (𝑃)
• 𝑃 is a vector that contains all the unknown parameters related to the BTC model • ‖⋅‖ is the Euclidean norm or 2-norm.

• 𝑗 𝑒𝑥𝑝 (𝑡) are the experimental current density.

• 𝑗 𝑠𝑖𝑚 (𝑡, 𝑃) are the simulated current density computed using the bipolar charge transport model in terms of P.

• 𝐶 1 (𝑃) represents the relative error between 𝑗 𝑒𝑥𝑝 and 𝑗 𝑠𝑖𝑚 (𝑃), so then 𝐶 1 (𝑃) is a unit-less scalar.

• 𝑛 𝑒𝑥𝑝 are the experimental data charge density.

• 𝑛 𝑠𝑖𝑚 (𝑃) are the simulated charge density computed using the bipolar charge transport model in terms of parameters P.

• 𝐶 2 (𝑃) represents the relative error between 𝑛 𝑒𝑥𝑝 and 𝑛 𝑠𝑖𝑚 (𝑃), so then 𝐶 2 (𝑃) is a unit-less scalar.

Chapter IV

123

• 𝑙𝑏 and 𝑢𝑏 are lower and upper bounds of the parameters respectively.

IV.5 Parameters optimization using TRRA IV.5.1 Influence of the cost function on the optimization outputs

In the following section, the three cases will be examined and compared to choose the most helpful cost function that can provide the best optimal solution compared to other cost functions.

IV.5.1.1 Cost function using current density

In this section, TRRA is used to calculate the unknown parameters to achieve the minimum squared difference between the experimental and simulated current density using three different applied fields (-20, -40 and -60 𝑘𝑉/𝑚𝑚). Consequently, the cost function to be minimized is defined as: When 𝐶 1 (𝑃) was minimized, TRRA produced an optimal set of parameters defined as 𝑃 1 . Figure 44 shows that 𝑃 1 was capable of producing a very good match between experimental and simulated current density for the three fields under consideration. The simulated current does not reach a stationary state for all considered fields, which is consistent with the experimental measurements. According to Figure 45, if the charge density is simulated using 𝑃 1 (parameters obtained by optimizing only current density data (𝐶 1 )), a significant difference could be noticed between experimental and simulated outcomes. Thus, 𝑃 1 is found as a local optimal solution that only fits the experimental current measurements and contradicts the space charge experimental measurements. However, our target is to find a global or unique set of parameters that is able to fit both current and charge measurements at the same time. It is worth mentioning that charge density data may provide more information than current measurements as the space charge measurements are offered with respect to time and position, whereas the current measurements are only with respect to time. For this reason, the next section will focus on optimizing the parameters by minimizing the difference between experimental and simulated charge density instead of current density.

𝐶 1 (𝑃) = ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -20𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -20𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -60𝑘𝑉⋅𝑚𝑚 -1

IV.5.1.2 Cost function using charge density

Here TRRA is used to minimize the squared difference between the experimental and simulated charge density using three different electric fields. Thus, the cost function to be optimized is defined as:

𝐶 2 (𝑃) = ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -20𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -20𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -60𝑘𝑉⋅𝑚𝑚 -1
Chapter IV 125

The optimal set of parameters found by optimizing 𝐶 2 (𝑃) is defined as 𝑃 2 and the obtained results are displayed in Figure 46. Based on this figure, there is a good correlation between experimental and simulated data for all applied fields, with the simulated patterns being compatible with the experimental ones. Furthermore, most of the properties highlighted in the experiment are reflected by the model, such as the fast injection of negative charges at short period and then the advantage of positive charges inside the bulk after 500 seconds. Same as what we have done previously, let us check the current density obtained by simulating BCT model using 𝑃 2 , which are represented in Figure 47.

For all considered applied fields, the difference between the experimental and simulated current density is found to be more than a decade. Despite the fact that charge density provides more information than current measurements, when the current density is simulated using 𝑃 2 , a significant mismatch could be noticed between experimental and simulated outcomes (Figure 47).

Thus, 𝑃 2 is found as a local optimal solution that only fits the experimental charge density measurements and contradicts the current experimental measurements.

Chapter

IV.5.1.3 Optimization using both current and charge densities

As mentioned previously, our main target is to find a unique set of parameters able to fit both current and charge experimental data simultaneously. However, according to the results represented in case 1 and case 2, it has been proved that 𝑃 1 and 𝑃 2 can produce a good correlation between experimental and simulated current and charge densities individually. On the other side, when 𝑃 1 is used to simulate the charge density, or when 𝑃 2 is used to simulate current density, an apparent discrepancy between experimental and simulated results could be observed. For this reason, TRRA is used in this section to minimize the squared difference between the experimental and simulated of both charge and current densities. Both cost functions will be combined in a single cost function 𝐶 3 (𝑃). In order to combine 𝐶 1 and 𝐶 2 in a single cost function 𝐶 3 , both functions were normalized as you can see in Equations 35 and 36 (Unit-less outputs). The cost function to be optimized is defined as:

𝐶 3 (𝑃) = 𝐶 1 (𝑃) + 𝐶 2 (𝑃) 𝐶 3 (𝑃) = ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -20𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -20𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑗 𝑒𝑥𝑝 -𝑗 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑗 𝑒𝑥𝑝 ‖ -60𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -20𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -20𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 𝑒𝑥𝑝 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 𝑒𝑥𝑝 ‖ -60𝑘𝑉⋅𝑚𝑚 -1
The optimized set of parameters 𝑃 3 produced by minimizing 𝐶 3 (𝑃) are displayed in Table 12, 13 and 14. The comparison between the experimental measurements and the simulated results using Figure 48 shows that 𝑃 3 can produce a good match between experimental and simulated current density, for all considered fields. The difference in current density between the experimental and simulated results is less than a quarter of a decade. It is worth noting that it is extremely difficult to find a single set of parameters that can match numerous current measurements simultaneously, because the experimental current measurements are not stable, and the curve shape may change if the same experiment is repeated twice. This is one of the factors that influenced the correlation between experimental and simulated current. Chapter

IV.5.2 Parameter analysis

Section 5 shows that the optimized parameters 𝑃 3 are more reliable than 𝑃 1 and 𝑃 2 since 𝑃 3 can suit both current and charge densities. Thus, the discussion part will be focused on analyzing and interpreting the results provided by 𝑃 3 (Table 12, 13 and 14).

IV.5.2.1 Recombination rate

The recombination reduction pre-factors 𝑅 𝑖 is introduced in order to capture the effect of reduced Langevin recombination. Based on Table 12, the optimal values of the recombination reduction pre-factors (𝑅 1 = 𝑅 2 = 𝑅 3 = 1), prove that Langevin's recombination form fully describes the recombination process with no reduction is taking place. The recombination rate could be evaluated by substituting 𝑅 𝑖=1,…,3 in Equations 38, 39 and 40.

For recombination rate between mobile electrons (𝑒𝜇) and trapped holes (ℎ𝑡):

Chapter IV 129 𝑆 𝑒𝜇-ℎ𝑡 = 𝑅 1 𝜇 𝑒 𝜀 0 𝜀 𝑟 [START_REF] O'dwyer | The theory of electrical conduction and breakdown in solid dielectrics[END_REF] For recombination between mobile holes (ℎ𝜇) and trapped electrons (𝑒𝑡):

𝑆 𝑒𝑡-ℎ𝜇 = 𝑅 2 𝜇 ℎ 𝜀 0 𝜀 𝑟 [START_REF] Hamann | Electrical conduction mechanisms in solids[END_REF] For recombination between mobile electrons (𝑒𝜇) and mobile holes (ℎ𝜇):

𝑆 𝑒𝜇-ℎ𝜇 = 𝑅 3 (𝜇 𝑒 + 𝜇 ℎ ) 𝜀 0 𝜀 𝑟 [START_REF]Handbook of nanoceramics and their based nanodevices[END_REF] In our model, the recombination rate is approximately stable with respect to time and position.

Thus, the recombination rate is evaluated by computing its average value with respect to time and position. and trapped holes (black). In addition to the discussed field dependencies of the recombination rates, their temperature dependency also exists and follows the consideration of the mobility equation.

IV.5.2.2 Trapping and detrapping rates

Figure 51 represents the optimal values of trapping and detrapping rates at fields ranging from -20 𝑘𝑉/𝑚𝑚 to -60 𝑘𝑉/𝑚𝑚 found by TRRA. Based on Figure 51, both trapping and detrapping rates are found to be increasing with the increase of the electric field. Besides, the trapping rates of electrons and holes are greater than the detrapping rates at any electric field. Also, the trapping rate of electrons is much higher than that of holes. Moreover, the rise of the trapping and detrapping rates of electrons with respect to the electric field is insignificant compared to the increase of the rates for holes. . Where 𝑥 is the difference in y-coordinates for the two points, 𝑡 is the difference in x-coordinates for these two points. Since only few amount of injected charges are detected at the beginning of the experiment, the range of the color map is reduced in this part to [-1,1] 𝐶/𝑚 3 to have a pure representation of the injection process of negative and positive charges. Only the experiment with -20 𝑘𝑉/𝑚𝑚 applied electric field could be used in this case to estimate the mobility of carriers.

For the PEA method, the electric field is applied with a ramp of 0.25 𝑘𝑉/𝑠 and thus it needs around 10 𝑠 to achieve -20 𝑘𝑉/𝑚𝑚 applied electric field. Moreover, PEA needs more than 30 𝑠 to reach -40 or -60 𝑘𝑉/𝑚𝑚 applied electric field. However, as mentioned previously, we should use the data provided at a short period (𝑡 < 15 𝑠) in order to get a clear representation of carrier mobility.

For this reason, the experiment with -20 𝑘𝑉/𝑚𝑚 electric field is the only one that could contribute to estimating the experimental mobility of both carriers. It is worth mentioning that this kind of approximation is not totally accurate, but it could provide additional useful information, such as the range of mobility values for both carriers. Where 𝐸 0 is the initial applied electric field, and 𝜇 𝑒,ℎ 𝑒𝑥𝑝 is the approximated experimental mobility of electrons and holes. The results are displayed in Table 15.

On the other hand, the simulated mobility of electrons and holes obtained by the BCT model using

𝑃
Hence, experimental and simulated mobility values of electrons and holes are displayed in Table 15. According to Table 15, (mobility values obtained from the experimental measurements), the mobility of electrons is found greater than the mobility of holes for any applied field. These results are in agreement with the simulated mobility values that are displayed in Table 15 (mobility values obtained by simulating the BCT model using 𝑃 3 ), where the mobility of electrons is greater than that of holes. Besides, the simulated mobility of electrons is found to be sensitive with respect to electric field compared to the mobility of holes. For example, the mobility of electrons increased from 3.5× 10 -13 to 9× 10 -13 𝑚 2 . 𝑉 -1 . 𝑠 -1 , however, the mobility of holes increased slightly from 1.35× 10 -14 to 1.8× 10 -14 𝑚 2 . 𝑉 -1 . 𝑠 -1 . Furthermore, the experimental and simulated mobilities of electrons are found to be approximately the same.

IV.5.2.4 Initial charge conditions

The optimal values of the initial net charge density parameters are displayed in Table 14. While performing the space charge measurements (i.e., Experiment 1, Figure 40), the initial charge density inside the bulk is found to be increasing after performing a new experiment on the same LDPE sample (i.e., 𝛼 𝐸𝑥𝑝1 -60 < 𝛼 𝐸𝑥𝑝1 -40 < 𝛼 𝐸𝑥𝑝1 -20 -20 ). This may be due to the fact that two different discharging techniques were used for each experiment. For the charge space measurements, the sample was manually discharged by placing the sample between two aluminum discs. Whereas for the current measurements, the sample was discharged by grounding both electrodes without touching the sample.

IV.6 Conclusion

In order to improve the BCT model, TRRA has been applied to identify a unique optimal set of parameters that best match both current and charge density measurements at varying electric fields (-20, -40, -60 𝑘𝑉/𝑚𝑚). Three distinct cases were examined in this chapter, where each example assumed a specific cost function to highlight the difficulties of finding a unique parameter. First, the cost function considered only the current density data for the optimization process. The outputs of TRRA revealed a strong correlation between experimental and simulated current density. Furthermore, if the charge density is simulated using these parameters, a significant mismatch is noticed between the experimental and simulated charge density. In the second case, the cost function examined only the charge density data in the optimization process.

The outputs of TRRA provided a good match between experimental and simulated charge density, but does not fit the experimental current density. Finally, the third cost function was created to incorporate both charge and current data in the optimization process. The findings obtained by applying the third cost function demonstrated a strong correlation between experimental and simulated outputs for both current and charge density data. Then, all the obtained optimal parameters were analyzed and discussed. The optimal values of the recombination reduction prefactors prove that the recombination process is entirely characterized by Langevin's recombination form with no reduction taking place. The optimal values of the mobility indicate that the mobility of electrons is higher than that of holes, which is consistent with the literature. Furthermore, the The optimal values of trapping and detrapping rates revealed that both trapping and detrapping rates increase with the increase of the electric field. Besides, the trapping rate is found to be greater than the detrapping rate for both carriers at any field. Moreover, the trapping rate of electrons was found to be much greater than that of holes. Whereas the detrapping rates of electrons and holes are found to be approximately equal.

Finally, the initial charges inside the bulk were observed to be rising after performing a new experiment on the same LDPE sample, even though the sample was discharged for 24 hours. The results also anticipated that the discharge strategy influenced the discharging process. The manual discharging technique used for charge density measurements seems to be more effective than maintaining the sample between electrodes.

The obtained optimal parameters properly suited both current and charge densities at varying electric fields. As noted earlier, the major objective of this work is to find a unique set of parameters that match both current and charge using any experimental conditions. In this study, only experiments with varying electric fields were examined, however, if the polarity of the electric field is reversed, the parameters may not suit the experimental outputs anymore. Thus, our technique might be further enhanced by introducing additional experiments with alternative experimental conditions (i.e., opposite polarity or varying temperatures) or adding another source of measurements (i.e., electroluminescence).

V.1 Introduction

In the previous chapter, TRRA was used to find a unique set of parameters to fit both current and charge density data with varying electric fields. The obtained parameters produced a good correlation between experimental and simulated results. Unfortunately, the obtained parameters are not sufficient to predict the charge behavior for all experimental conditions. For this reason, the physical phenomena should be critically analyzed individually to understand their impact on the charge distribution. In this chapter, we will only focus on studying the charge injection phenomena taking place at the electrodes using optimization tools.

In practice, the effects of electrode materials on charge injection processes are often addressed and have become a growing area of interest due to their influence on the space charge behavior [1].

Early experimental studies [2], [3] indicated that the electrode materials have a significant effect on the magnitude of the conduction current and space charge distribution. The major difficulty arises from the fact that there exist only a little knowledge regarding the charge injection at the electrodes and the trapping processes taking place in the bulk of LDPE.

In the present study, the effect of four different electrode materials (i.e., gold (Au), aluminum (Al), semiconductor (Sc), and copper (Cu)), on the charge injection process will be investigated using the PEA technique. In order to reduce the influence of impurities such as antioxidants and cross-linking byproducts on space charge formation [3], additive-free LDPE was selected in the present study so that the effect of electrode materials on charge formation can be examined.

It is important to note that the work presented in this chapter is a preliminary study. It is the first step to better understand the injection mechanism at the surface electrode in order to enhance its numerical modeling behavior.

This chapter is divided into four sections. After this brief introduction, the second section will give a detailed explanation of the material preparation and the experimental protocol done in our laboratory. Besides, we will present the experimental measurements obtained by PEA using different surface electrodes. The third part will be concerned with optimizing the parameters using the experimental data, and finally, we will discuss and analyze the results.

Chapter V 140

V.2 Experimental measurements V.2.1 Material preparation and Experimental Protocol

In this chapter, four LDPE samples with the same characteristics are manufactured to perform the experimental measurements using different electrodes. Indeed, the LDPE films are prepared using the same technique mentioned in chapter 4. Each sample is metallized with different kinds of materials on the surface. All samples are 150 µm thick and have a diameter of 80 mm. Three different metallic electrodes (Au, Al and Cu) with a diameter of 10 mm and a 50 nm thickness were sputter-coated onto the LDPE samples. In the case of Sc, a film of 75 μm thickness with a diameter of 10 mm was manufactured in our laboratory using a heating press. The Sc was attached to the LDPE sample by adding some silicon oil between them to guarantee good contact.

The same experimental protocol is applied to the four different samples discussed previously The sample has been discharged 24 hours after each measurement to eliminate the charges in the bulk as much as possible. The samples have been placed accurately in the PEA cell to ensure that all measurements are carried out at the same location while applying different electric fields.

V.2.2 Charge density profiles for all experiments

After performing all the measurements displayed in the protocol, the experimental results are displayed in Figure 54, 3, 4 and 5. All of the cartography in this section depicts the variation of Chapter V 141 space charge behavior with respect to time and position under a DC stress at -40 and -60 𝑘𝑉/𝑚𝑚 (-6 𝑘𝑉 and -9 𝑘𝑉, respectively).

In the following we try to study how different types of surface electrodes affect the injection process. However, this study is difficult to be done based on experimental data provided by PEA method, which are partially linked to the limited spatial and temporal resolution. For this reason, we intend to compare the nature of different electrodes using only the charge distribution cartography. These cartographies provide a tendency of charge build-up at short time that could be used to analyze the impact of surface electrodes on the injection process. For -40 𝑘𝑉/𝑚𝑚 of applied field (left side of Figure 54), positive (red) and negative (blue) charges are detected immediately after the voltage is applied. The electrons could be observed with a strong and rapid injection directly after applying the voltage. Moreover, the mobility of electrons seems to be higher than that of holes. For a short time (t < 300 s), the number of negative charges is greater than the number of positive ones inside the bulk of LDPE sample. For a time greater than 600 s, positive charges predominate inside the bulk of the insulator, compared to negative charges.

V.2.2.1 Space charge distribution with gold electrodes (Au-LDPE-Au)

Furthermore, after 1500 s of polarization, the insulator seems to reach an equilibrium state in terms of the space charge distribution, since the net density of charges remains stable until the end of the experiment. For -60 𝑘𝑉/𝑚𝑚 of applied field (right side of Figure 54), the same behavior takes place inside the bulk as it does under -40 𝑘𝑉/𝑚𝑚 of applied field. Furthermore, the number of injected charges increases as the applied voltage is increased.

Chapter V 142 the mobility of electrons seems to be higher than that of holes. After 1000 seconds, the behavior of the insulator seems to reach an equilibrium state in terms of the space charge distribution and amount of injected charge (i.e., the number of positive and negative charges remains stable until the end of the experiment). Moreover, the number of positive charges in the bulk of LDPE with aluminum electrodes appears to be less than the positive charges in the bulk of LDPE with gold electrodes. For -60 𝑘𝑉/𝑚𝑚 applied field (right), the same behavior takes place inside the bulk as it does under -40 𝑘𝑉/𝑚𝑚 of applied field. Moreover, the number of injected electrons significantly increases with the increase in the applied field.

V.2.2.2 Space charge distribution with Aluminum electrodes (Al-LDPE-Al)

V.2.2.3 Space charge distribution with copper electrodes (Cu-LDPE-Cu)

According to Figure 56, for -40 𝑘𝑉/𝑚𝑚 (Left) applied field, the negative charge could be seen with a strong and rapid injection directly after applying the voltage. Indeed, the behavior of negative charges is compatible with the behavior of LDPE with gold and aluminum electrodes at short time. Moreover, the mobility of electrons seems to be higher than that of holes. At short time (𝑡 < 300 seconds), the number of negative charges is greater than the number of positive ones inside the bulk of LDPE sample. Furthermore, for a time greater than 500 𝑠, positive charges have the advantage within the bulk of the insulator, compared to negative charges. The behavior of positive charges after 500 s is identical to the positive charges behavior for LDPE with gold electrodes. Furthermore, for 5400 s of polarization, the insulator did not reach an equilibrium state in terms of the space charge distribution, where the net density of charges drops with respect to time. The equilibrium state was not reached maybe due to the impact of the trapped and detrapping phenomena or another microscopic phenomena taking place inside the bulk such as recombination between carriers. For -60 𝑘𝑉/𝑚𝑚 applied field (right), the same behavior takes place inside the bulk as it does under -40 𝑘𝑉/𝑚𝑚 of applied field. Moreover, the number of negative and positive charges significantly increases with the increase in the applied field. According to Figure 57, for -40 𝑘𝑉/𝑚𝑚 (Left) applied field, the negative charge could be seen with a strong and rapid injection at the beginning of the experiment. Moreover, the mobility of Chapter V 144 electrons seems to be higher than that of holes. For short time (t < 300 seconds), the number of negative charges is greater than the number of positive charges. For a time greater than 600 s, positive charges have the advantage inside the bulk of the insulator compared to negative charges.

V.2.2.4 Space charge distribution with Semiconductor electrodes (Sc-LDPE-Sc)

Furthermore, the insulator did not reach an equilibrium state in terms of the space charge distribution as the net density of charges drops after time passes.

For -60 𝑘𝑉/𝑚𝑚 applied field (right), the same behavior takes place inside the bulk as it does under -40 𝑘𝑉/𝑚𝑚 of applied field. Moreover, the number of negative and positive charges significantly increases with the increase in the applied field.

V.2.3 Qualitative analysis

Experimentally, the effect of surface electrodes on the injection process could be observed at short time after applying the electric field. After long time of polarization, it is difficult to understand the effect of surface electrodes on the injection process because many additional phenomena took place inside dielectrics, such as recombination, trapping and detrapping phenomena.

The measurements demonstrated that for both electrons and holes, semi-conductor injects the maximum amount of charges compared to other electrodes, then the copper was found in the second stage. Furthermore, the gold and aluminum surfaces were found with the minimum charge injection compared to semi-conductor and copper. Besides, for gold and aluminum measurements, it is difficult to identify which material inject more than the other using PEA method. Thus, from the carrier injection point of view, the experimental measurements indicated that for electrons and holes, the increasing order depending on material type is as follows: 𝐴𝑙 ≈ 𝐴𝑢 < 𝐶𝑢 < 𝑆𝑐.

V.3 Optimization using all experimental data

After presenting the experimental measurements and analyzing the obtained results, now TRRA will be applied to estimate the injection barrier heights of each surface electrode, as well as all the other unknown parameters related to BCT model. The aim of this part is to analyze the capability of electrode materials to inject charges into dielectrics by comparing two different approaches, experimentally and simulation.

V.3.1 Unknown parameters

The units, symbols, and bounds of all the considered unknown parameters are displayed in As mentioned previously, the goal of this section is to understand the effect of the surface electrode on the injection process. This could be achieved by computing the values of the barrier height of injection of each of the surfaces used in this section (Au, Al, Cu, and Sc). Thus, four additional unknown parameters are introduced to the optimization process. The unknown injection barrier height parameters are represented in Table 16. According to the experimental protocol (Figure 53) used in this chapter, two different electric fields (-40 and -60 𝑘𝑉/𝑚𝑚) are applied to the same LDPE sample with different surface electrodes (Au, Al, Cu, and Sc). As in Chapter 4, the trapping and detrapping parameters are supposed to vary with the variation of the electric field. For this reason, the trapping and detrapping parameters to be found are as follows:

Unit 𝑠 -1 𝑠 -1 Range [0 -1] [0 -1] Chapter V 146
The initial charges inside the bulk are represented by the parameters 𝛼 𝑖 and 𝛽 𝑖 (see Figure 58), where 𝑖 = 𝐴𝑢, 𝐴𝑙, 𝐶𝑢 or 𝑆𝑐. Based on the results obtained in Chapter 4 (Table 14), the initial charges in virgin LDPE samples were found to be approximately equal to zero, thus, 𝛽 𝑖 = 0 𝐶/𝑚 3 for any surface electrode. Therefore, 4 parameters are introduced to the optimization process to express the initial charges before starting the experiments of -40 𝑘𝑉/𝑚𝑚 electric field, that are, 𝛼 𝐴𝑢 , 𝛼 𝐴𝑙 , 𝛼 𝐶𝑢 and 𝛼 𝑆𝑐 (Figure 58 and Table 19). In total, 27 unknown model parameters are included in the optimization process. The starting point of the optimization process are the optimal parameters obtained in Chapter 4. This was done to enhance the robustness and the efficiency of the optimization algorithm.

In this chapter, TRRA is used to find the unknown parameters that minimize the cost function 𝐶(𝑃) combining all the experiments.

The cost function to be optimized is defined as:

Chapter V 147 𝐶(𝑃) = ‖𝑛 exp 𝐴𝑢 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 exp 𝐴𝑢 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 exp 𝐴𝑢 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 exp 𝐴𝑢 ‖ -60𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 exp 𝐴𝑙 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 exp 𝐴𝑙 ‖ -60𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 exp 𝐴𝑙 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 exp 𝐴𝑙 ‖ -60𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 exp 𝐶𝑢 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 exp 𝐶𝑢 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 exp 𝐶𝑢 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 exp 𝐶𝑢 ‖ -60𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 exp 𝑆𝑐 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -40𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 exp 𝑆𝑐 ‖ -40𝑘𝑉⋅𝑚𝑚 -1 + ‖𝑛 exp 𝑆𝑐 -𝑛 𝑠𝑖𝑚 (𝑃)‖ -60𝑘𝑉⋅𝑚𝑚 -1 ‖𝑛 exp 𝑆𝑐 ‖ -60𝑘𝑉⋅𝑚𝑚 -1 V.3.

Charge distribution: experimental vs simulation

In this section, the experimental results obtained by PEA method are presented with the simulated results obtained by BCT model using the optimal parameters obtained by TRRA (Figure 59, 60, 61 and 62). All the obtained cartographies and parameters will be compared and discussed in the following sections. electrodes. The model was able to reproduce most of the characteristics presented in the experiment, such as the dominance of the negative charges, due to strong injection at short time (𝑡 < 300 𝑠). The advantage of positive charges for time greater than 600 s inside the bulk of the LDPE sample, compared to negative ones could also be observed. Moreover, the mobility of negative charges seems to be greater than that of positive ones at short time.

For gold electrodes, both experimental and simulated outputs show that insulator seems to reach an equilibrium state in terms of the space charge distribution for 𝑡 > 1000 𝑠.

For LDPE with aluminum electrodes, the model was able to reproduce most of the characteristics presented in the experiment, such as, the small quantity of positive charges injected compared to the quantity of negative charges. Also, the advantage of negative inside the bulk of insulation Chapter V 150 compared to positive charges at any time of polarization. Furthermore, for both experimental and simulated outputs, the insulator reached an equilibrium state in terms of the space charge distribution where the net density of charges keeps stable for 𝑡 > 1000 𝑠.

For LDPE samples with copper and semiconductor electrodes, the experimental measurements show that the insulator did not reach an equilibrium state in terms of the space charge distribution where the net density of charges keeps dropping inside the bulk as the time passes. On the other hand, the BCT model was not able to reproduce this decay of the net density of charges (i.e., after 1000 s, the simulated net density of charges remains stable).

V.3.3.2 Optimal parameters

The optimized set of parameters 𝑃 ⋆ produced by minimizing 𝐶(𝑃) are displayed in Table 20, 21, 22 and 23. All the optimal parameters are presented but the discussion will focus on the parameters representing the barrier height of injection and the initial net density of charges. These parameters are the ones linked to the nature of surface electrode. Based to the results displayed in Table 20, the barrier height for injecting positive charges is found to be greater than the barrier height for injecting negative charge for all the examined surface electrodes. According to the modified Schottky law (i.e., Equation 1), it is known that if the injection barrier height parameter increases, the quantity of injected carriers decreases, and the inverse is also true. Therefore, the quantity of injected electrons is found to be more than that of the holes through all electrodes.

Moreover, these optimal parameters can be used to compare the quantity of injected carriers through surface electrodes of different nature.

For negative charges:

The semi-conductor is found to be the material that injects the highest number of negative charges compared to all other surfaces electrodes (i.e., lowest barrier height compared to other materials, 𝑤 𝑒 𝑆𝑐 = 0.980 𝑒𝑉). Then the copper (𝑤 𝑒 𝐶𝑢 = 1.116 𝑒𝑉) was found in the second stage, where it injects more electrons than gold (𝑤 𝑒 𝐴𝑢 = 1.133 𝑒𝑉) and aluminum (𝑤 𝑒 𝐴𝑙 = 1.140 𝑒𝑉). Finally, the aluminum surface injects the minimum quantity of negative charges compared to all other surfaces. From the carrier injection point of view, the results indicate that for electrons, the increasing order depending on material type is as follows: 𝐴𝑙 < 𝐴𝑢 < 𝐶𝑢 < 𝑆𝑐.

For positive charges: Based on the optimal values of the injection barrier height, it is obvious that semi-conductor is the material that injects the highest number of positive charges compared to all other surfaces (i.e., lowest barrier height compared to other materials, 𝑤 ℎ 𝑆𝑐 = 1.000 𝑒𝑉). Then in the second stage, we find the copper (𝑤 ℎ 𝐶𝑢 = 1.132 𝑒𝑉), which injects more positive charges than gold and aluminum. Finally, the injection of positive charges through gold (𝑤 ℎ 𝐴𝑢 = 1.170 𝑒𝑉) is more than the injection through aluminum (𝑤 ℎ 𝐴𝑙 = 1.194 𝑒𝑉). Thus, from the carrier injection point of view, the increasing order depending on material type for holes is as follows: 𝐴𝑙 < 𝐴𝑢 < 𝐶𝑢 < 𝑆𝑐. These results are similar to the classification obtained for negative carriers.

Experimentally, from the carrier injection point of view, the qualitative analysis presented in subsection V.2.3 indicated that for both electrons and holes, the increasing order depending on material type is as follows: 𝐴𝑙 ≈ 𝐴𝑢 < 𝐶𝑢 < 𝑆𝑐. Therefore, the optimal parameters obtained by TRRA are found to be compatible with the experiment measurements.

Chapter V 152

V.3.3.2.2 Initial charge density

The optimal values of the initial net charge density parameters are presented in Table 21. The obtained results reveal that only a small quantity of initial charges is detected inside the bulk of the sample for LDPE with metal electrodes (i.e., 𝛼 𝐴𝑢 ≈ 𝛼 𝐴𝑢 ≈ 𝛼 𝐴𝑢 ≈ 0.1 𝐶/𝑚 3 ). The amount of initial charges detected in LDPE with semi-conductor electrodes (i.e., 𝛼 𝑆𝑐 = 4.7 𝐶/𝑚 3 ), is significantly greater than the charges detected in LDPE with metal electrodes. These results suggest that the discharging of a sample is easier when using metal electrodes instead of semi-conductor. Thanks to TRRA, we are now able to have a better knowledge of the influence of surface materials on the discharging phenomena. However, additional experiments should be performed to confirm the validity of these results.

V.4 Conclusion

TRRA has been applied to identify the optimal set of parameters that best fit charge density measurements with four different surface electrodes at varying electric fields. The optimal parameters provided a good correlation between experimental and simulated charge density data for -40 𝑘𝑉/𝑚𝑚 and -60 𝑘𝑉/𝑚𝑚 electric fields, for all considered surface electrodes. Besides, the model was able to reproduce most of the characteristics presented in the experimental measurements. The optimal values of the injection barrier height parameters were found to be correlated with the experimental measurements, which indicate that for both electrons and holes, the order is as follows: 𝑆𝑐 > 𝐶𝑢 > 𝐴𝑢 > 𝐴𝑙. Furthermore, the parameters representing the injection barrier heights of 𝐴𝑢 and 𝐴𝑙 were found near to each other, which is similar to the results observed experimentally.

The optimal parameters representing the initial charges inside the bulk show that the discharging of a sample is easier when using metal electrodes instead of semi-conductor.

As mentioned previously, the obtained optimal parameters properly suited all experiments with different surface electrodes and with varying electric fields. Unfortunately, the obtained results do not yet match the parameters that we are intending to find. For example, the trapping, detrapping, and injection parameters were found to be different from the values obtained in Chapter 4. Thus, it is a challenging task to find a unique set of parameters that can fit several experiments with different experimental conditions (e.g., different surface electrodes, varying electric fields, or varying temperatures). This is one of the weaknesses of the BCT model that may be enhanced in our future work. rate, de-trapping rate, and recombination coefficients. Most of these parameters cannot be predicted by independent experiments. For this reason, the numerical modeling approach was coupled with the experimental approach in this study, which provided a lot of beneficial information that contributed to estimating the unknown parameters.

To enhance the BCT model and to be able to estimate the behavior of dielectrics under DC stress using any experimental conditions, one must find a unique set of parameters that can always provide a good correlation between experimental and simulated data. Therefore, the principal objective of this study was to provide an original approach that contributes to identifying this unique set of parameters. Consequently, optimization tools have been introduced to find the best set of parameters that minimize the sum of the squares of the deviations between experimental and simulated data. and weak points of each algorithm. In order to select the most reliable algorithm, a trivial example was employed to examine all the considered algorithms. TRRA was found to be the most suitable algorithm for this type of problem. Compared to the other methods, TRRA had the best convergence rate and the shortest computation time.

TRRA has been applied to identify the optimal set of parameters that best match both current and charge density measurements at varying electric fields (-20, -40, -60 kV/mm). Three distinct cases were examined, where each example assumed a specific cost function. This chapter highlights the difficulty of finding a unique set of parameters that can suit several experiments. For example, if the cost function considered only the current density data for the optimization process. The outputs of TRRA revealed a strong correlation between experimental and simulated current density.

Furthermore, if the charge density was simulated using the obtained parameters, a significant mismatch was noticed between the experimental and simulated charge densities. Thus, a multiobjective function was formed, aiming at incorporating both charge and current data in a single cost function. The findings obtained by the new cost function demonstrated a strong correlation between experimental and simulated outputs for both current and charge density data at varying electric fields. Nevertheless, only experiments with varying electric fields were examined in this study; however, if the polarity of the electric field is reversed or if the temperature changed, the parameters may not suit the experimental outputs anymore. Hence, our technique might be further enhanced by introducing additional experiments with alternative experimental conditions (i.e., opposite polarity or varying temperatures) or by adding another source of measurements (i.e., electroluminescence).

In addition, the physical phenomena in the BCT model were analyzed based on the obtained optimal parameters. For example:

• The optimal values of the recombination reduction pre-factors proposed that the recombination processes are entirely characterized by Langevin's recombination form with no reduction taking place.

• The optimal values of the mobility indicate that the mobility of electrons is higher than that of holes, which is consistent with the literature. Moreover, the experimental and simulated mobility values were found quite near to each other and varied in the same range.

Conclusion

159

• The optimal values of trapping and detrapping rates revealed that both trapping and detrapping rates increase with the increase of the electric field. Besides, the trapping rate was found to be greater than the detrapping rate for both carriers at any given field.

• Finally, the initial charges inside the bulk were observed to be rising after performing a new experiment on the same LDPE sample, even though the sample was discharged for 24 hours.

The results also anticipated that the technique used for discharging the sample might influence the discharging process.

All the physical phenomena were analyzed based on the simulation approach obtained by the BCT model, considering the model weaknesses that may affect the reliability of the outcomes (i.e., some features are neglected in the model for simplicity reasons).

Finally, TRRA has been applied to better understand the influence of surface electrodes on the injection barrier height parameters using four different surface electrodes (i.e., 𝐴𝑢, 𝐴𝑙, 𝐶𝑢, and 𝑆𝑐) at varying electric fields (i.e., -40 𝑘𝑉/𝑚𝑚 and -60 𝑘𝑉/𝑚𝑚). The obtained optimal parameters provided a good correlation between experimental and simulated charge density data for all considered surface electrodes.

Thanks to optimization techniques, the optimal values of the injection barrier height parameters were found to be compatible with the experimental measurements, which indicate that for both electrons and holes, the order is as follows: 𝑆𝑐 > 𝐶𝑢 > 𝐴𝑢 > 𝐴𝑙.

Moreover, the mobility of electrons was found to be greater than that of holes, and the trapping and detrapping rates increased with the increase in electric field. Hence, the optimal parameters of the mobility, trapping, and detrapping were found to be consistent with the results obtained in Chapter 4.

Finally, the optimal parameters representing the initial charges inside the bulk show that the discharging of a sample is easier when using metal electrodes instead of semi-conductor.

The obtained optimal parameters properly suited all experiments with different surface electrodes and with varying electric fields. Unfortunately, the obtained results do not yet match the parameters that we intend to find. For example, the values of trapping, detrapping, and injection parameters were found to be different from the values obtained in Chapter 4. Thus, it is challenging to find a unique set of parameters that can fit several experiments with different experimental conditions (e.g., different surface electrodes, varying electric fields, or varying temperatures, opposite polarity). This is one of the weaknesses of the BCT model that may be enhanced in our future work. Many additional features could be added to the BCT model that could contribute to optimizing it by finding a unique set of parameters that suit any experiment (i.e., diffusion, depolarization, etc..). Additional experiments may also contribute Conclusion

160
to construct an equation able to characterize the trapping and detrapping rates in terms of the electric field.

Indeed, this work is only a key step in the development of a model capable of predicting the charge distribution in a solid dielectric. The following are some of the future perspectives for this work:

• In any modeling process, sensitivity analysis approaches are critical for identifying the influence of the inputs on the outcome. This assists in the identification of the most important experimental conditions to consider during the optimization process. One of the most efficient global sensitivity analysis methods that could be used in this study is the Sobol Sensitivity Analysis method. Sobol sensitivity analysis can be used to reveal the most essential experiments to be evaluated in the optimization process, rather than employing a large number of experiments to find a unique set of parameters.

• Introducing another source of experiments into the optimization process, such as electroluminescence measurements. This can provide additional information that could help in finding the most accurate set of parameters that can fit any experimental measurement, especially recombination rates.

• In this study, we considered several experiments with varying electric fields. In future work, adding experiments with varying temperatures or reversing the electric field polarity can be helpful and may contribute to finding a unique set of parameters.

• In this study, the depolarization (after removing the electric field) process was not considered in the optimization process. Indeed, the data from the depolarization period could provide additional important information that could enhance the behavior of the optimization algorithm toward the global optimal solution.

• The optimization technique proposed in this study can be further improved by using meta-model approaches (e.g., kriging or SVM). These techniques could make it possible to replace the computationally expensive cost function by a much less expensive one (surrogate modeling).

These models will allow further sensitivity studies and global minimization approaches (e.g., genetic algorithm) relying on a large number of function evaluations to be performed.

Figure 1 :

 1 Figure 1: Some defects encountered in low-density polyethylene. a) conformational defects in the amorphous zone,b) examples of unsaturated species at the origin of the chemical disorder.[24]

Figure 2 :

 2 Figure 2: Simplified description of the polyethylene band gap

Figure 3 :

 3 Figure 3: Ethylene molecule, monomer of polyethylene.

Figure 4 :

 4 Figure 4: Chemical structure of polyethylene

Figure 5 :

 5 Figure 5: Schematic energy band diagram of Schottky emission in metal-insulator structure.

Figure 6 :

 6 Figure 6: Schematic energy band diagram of Fowler-Nordheim emission in metal-insulator structure I.3.1.3 Thermionic-field emission Thermionic-field emission occurs between field emission and Schottky emission. In this condition, the tunneling electrons should have an energy level between the Fermi level of metal and the conduction band edge of the dielectric. The schematic energy band diagram of thermionic-field

Figure 7 :

 7 Figure 7: Schematic energy band diagram of thermionic-field emission in metal-insulator structure.• To summarize what has been stated so far, charges can be produced in dielectric films by three different conduction mechanisms: ionization, irradiation, or injection at electrodes. In our case, the charges are assumed to be produced only by injection at both electrodes,

Figure 8 :

 8 Figure 8: Comparison of thermionic-field emission, thermionic emission, and field emission I.3.2 Charge transport mechanism The charge transport mechanisms depend on the electrical properties of the dielectric itself. The trap energy level in the dielectric films is the most important parameter in this type of conduction process. The charge transport mechanisms include: • Poole-Frenkel emission • Hopping conduction • Ohmic conduction • Space charge-limited conduction • Diffusion conduction I.3.2.1 Poole-Frenkel Emission.

Figure 9 :

 9 Figure 9: Schematic energy band diagram of Poole-Frenkel emission in metal-insulator structure.

Figure 10 :

 10 Figure 10: Energy band diagram of hopping conduction in metal-insulator structure.

Figure 11 :

 11 Figure 11: Energy band diagram of ohmic conduction in metal-insulator structure.

Figure 12 :

 12 Figure 12: Schematic representation of the conduction and trapping/de-trapping model for LDPE. 𝑤 𝑒 𝑎𝑛𝑑 𝑤 ℎ are injection barriers for electrons and holes; 𝑅 0 , 𝑅 1 , 𝑅 2 , 𝑎𝑛𝑑 𝑅 3 are reduction pre-factors multiplied by the Langevin's recombination reduction pre-factor, 𝐵 𝑒 𝑎𝑛𝑑 𝐵 ℎ represents the trapping parameters, 𝐷 𝑒 𝑎𝑛𝑑 𝐷 ℎ represents the detrapping parameters; 𝑎 𝑒 𝑎𝑛𝑑 𝑎 ℎ represents the distance between two shallow traps, 𝑤 𝑚𝑜 𝑒 and 𝑤 𝑚𝑜 ℎ are electron and hole depth for a single shallow trapping level.

  [START_REF] Chen | A new model for surface potential decay of corona-charged polymers[END_REF] proposed a modified bipolar charge injection model to account for surface charge decay of corona-charged polymeric materials. The model can take care of field dependent carrier mobility and readily explains the surface potential decay with different charging times. The results provided by corona-charged samples shows a double injection from both electrodes. The new model reveals that the surface potential cross-over phenomenon can occur under bipolar charge injection.F. Baudoin et al. (2011-2012):The authors suggest a charge transport model that considers the electroluminescence phenomena (EL) in polyethylene films under ac stress[START_REF] Baudoin | Modeling electroluminescence in insulating polymers under ac stress: effect of excitation waveform[END_REF]. Charge injection/extraction, transport, and recombination are also included in the model. The following are some of the modifications made by Baudoin:

Figure 13

 13 Figure 13 shows an example of a one-dimensional optimization problem.

Figure 13 :

 13 Figure 13: A one-dimensional optimization problem. Note that the minimum is the best in the feasible set, lower points may exist outside the feasible region.

  the function's highest and lowest values, either within a specific range (local) or over the entire domain (global). Indeed, we can simply demonstrate that the problems min 𝑥∈𝒳 𝑓(𝑥) and max 𝑥∈𝒳 -𝑓(𝑥) (with or without constraints) are equivalent since they contain the same set of solutions, such that:

•

  Inflection: Figure 14 also shows an inflection point where the derivative is zero, but the point does not locally minimize 𝑓. An inflection point is where the sign of the second derivative of 𝑓 changes, which corresponds to a local minimum or maximum of 𝑓′. • Global: Most optimization algorithms work by first trying to locate any feasible solution, and then attempting to find another (better) feasible solution that minimizes the value of the objective function. The global minimum of an objective function is achieved at the global solution, as presented in Figure 14. A point 𝑥 ⋆ is a global minimum if: 𝑓(𝑥 ⋆ ) ≤ 𝑓(𝑥), ∀𝑥 ∈ ℝ 𝑛 .

Figure 14 :

 14 Figure 14: Examples of critical points of interest to optimization algorithms (where the derivative is zero) on a univariate function

Figure 15 :

 15 Figure 15: Examples of the necessary but insufficient conditions for strong local minima.

  it may be used to describe all of the previous examples. In other words, LSO problems might be linear, nonlinear, constrained, or unconstrained. The LSO technique is a classic regression analysis strategy for approximating the solution of overdetermined systems (sets of equations with more equations than unknowns) by minimizing the sum of the squares of the residuals in the outcomes of each individual equation. Least-squares problems occur while fitting a parameterized mathematical model to a set of data points by minimizing an objective defined as the sum of the squares of the model function and data point errors. There are two types of least-squares problems: linear and nonlinear least-squares, based on whether the residuals are linear in all unknowns.

Algorithm 4 :

 4 Newton's methodInitial conditions: 𝑥 0 ∈ ℝ 𝑛 1: while stop criteria is not achieved, do 2: Calculate Newton's direction 𝑠 𝑘 the solution of 𝛻 2 𝑓(𝑥 𝑘 ) ⋅ 𝑠 𝑘 = -∇𝑓(𝑥 𝑘 )3:𝑥 𝑘+1 = 𝑥 𝑘 + 𝑠 𝑘 4: 𝑘 ← 𝑘 + 15: return 𝑥 𝑘Chapter II 69 on derivatives or approximate them via finite differences are insignificant. In such circumstances, the optimization problem is called derivative-free optimization.

  These algorithms are used to minimize the sum of the squares of the deviations between experimental data and simulation data. Two experimental sources could be employed in our case, namely charge density and current density measurements. Charge density is measured by the PEA method along with what is known as charging current measurement. On the other hand, a bipolar charge transport model developed for LDPE is used to generate the simulated data.

Figure 16 :

 16 Figure 16: Principle of the optimization technique. P is the set of unknown parameters, x represents the displacement, t represents the time, 𝑓 𝑒𝑥𝑝 are the experimental data, 𝑓 𝑠𝑖𝑚 are the simulated data, EC are the experimental conditions.

  optimization algorithms have been proposed for solving blackbox cost functions. Such algorithms do not expect first or second-order derivatives to be analytically accessible. The black box optimization algorithms are composed of two categories: the direct search algorithms and stochastic algorithms. In the following sections, 3 types of direct search algorithms will be addressed: Trust Region Reflective, Levenberg-Marquardt, and Nelder-Mead. Besides, 2 types of stochastic algorithms are considered: Genetic and Particle Swarm algorithms.

1 2𝑟

 1 𝑇 (𝑃)𝑟(𝑃), 𝛻𝐶(𝑃) = 𝐽 𝑟 (𝑃)𝑟(𝑃) 𝛻 2 𝐶(𝑃) = 𝐽 𝑟 (𝑃)𝐽(𝑃) + 𝐻(𝑃)𝑟(𝑃)

Figure 17 :

 17 Figure 17: Example of a simplex in ℝ 2 Construct a symmetric point 𝑃 𝑠 = 𝑃 𝑔 + (𝑃 𝑔 -𝑃 2 ) of 𝑃 2 with respect to the center of gravity 𝑃 𝑔 , as represented in Figure 18.

Figure 18 :

 18 Figure 18: Reflection operation of the simplex if 𝐶(𝑃 𝑠 ) < 𝐶(𝑃 0 ), then enlarge the simplex in the direction of 𝑃 𝑠 , such that, 𝑃 𝑎 = 𝑃 𝑔 + 2(𝑃 𝑔 -𝑃 2 ), as shown in Figure 19. if 𝐶(𝑃 𝑎 ) < 𝐶(𝑃 𝑠 ) then 𝑃 2 is eliminated and 𝑃 𝑎 becomes the new vertex of the new simplex

Figure 22 :

 22 Figure 22: The description of one iteration of the genetic algorithm's evolution which operates in four stages: Evaluation, where we calculate the fitness of each parameter randomly generated; Selection, where it chooses a relatively fit subset of individuals for breeding; Crossover, where it recombines pairs of breeders to create a new population; and Mutation, where it potentially modifies portions of new chromosomes to help maintain the overall genetic diversity. Arrows in the figure indicate the transitions into the next genetic operation within one generation.

Figure 23 :

 23 Figure 23: Updating Particle Position and Velocity

Where

  𝑛 𝑆𝑖𝑚 (P) are the simulated charge density by the BCT model in terms of unknown parameters P and 𝑛 𝐸𝑥𝑝 are the experimental charge density obtained by PEA.

Figure 24 shows

 24 Figure 24 shows the experimental space charge behavior as a function of time and position with 500 𝑠 of charging cycle on a 200 𝜇𝑚 thick LDPE film containing some antioxidants, with -30 𝑘𝑉/𝑚𝑚 applied field at 40 °𝐶. The x-axis represents the polarization time (hours) of the experiment, and the y-axis represents the position (μm). The color bar represents the quantities of electrons and holes in terms of 𝐶/𝑚 3 . The blue color represents the electrons and the red color represents the holes. In order to avoid the representation of the capacitive and image charges beside the electrodes, only the data in the range of [25 µm -175 µm] are presented.

Figure 24 :

 24 Figure 24:Experimental net density of charge. Applied field: -30 kV/mm, 500s of polarization, at 40 °𝐶 Sample thickness: 200 µm. 25 µm beside electrodes are neglected, color bar provides charge density scale in 𝐶 ⋅ 𝑚 -3

\

Figure 25 :

 25 Figure 25: Simulated Charge density obtained by BCT model. TRRA: simulated by using 𝑃 𝑇𝑅𝐴 , LM: simulated by using 𝑃 𝐿𝑀 , NM: simulated using 𝑃 𝑁𝑀 , GA: simulated by using 𝑃 𝐺𝐴 , and PS: simulated by using 𝑃 𝑃𝑆

Figure 26 :Figure 27 :Figure 28 :

 262728 Figure 26: Press equipment used for preparing LDPE films from pellets

Figure 29 :

 29 Figure 29: Sputtering machine used for metallization

Figure 31 :

 31 Figure 31: Current measurement cell

  represents the absolute value of the current density with respect to polarization time under a -60 𝑘𝑉/𝑚𝑚 (-9 𝑘𝑉) DC stress with 3 hours of charging cycle on LDPE film of thickness 𝐿 = 150 𝜇𝑚. The x-axis represents the polarization time (seconds) of the experiment, and the y-axis represents the current density (A/m 2 ). It is clear that the charge density does not reach a stable state where it keeps decreasing, even after 3 hours of polarization.

Figure 32 :

 32 Figure 32: Absolute experimental current density. Applied field: -60 𝑘𝑉/𝑚𝑚 at room temperature, 3 hours of polarization, LDPE coated with a gold layer at both sides. Sample thickness: 150 μm.

Figure 33 :

 33 Figure 33: Principle of PEA method

Figure 34 :Figure 35 :

 3435 Figure 34: PEA device

Figure 36 :

 36 Figure 36: Experimental protocol of polarization/short-circuit (150 µm thick LDPE film).

Figure 37 shows

 37 Figure 37 shows the variation of space charge behavior with electric field as a function of time and position at -20, -40, and -60 𝑘𝑉/𝑚𝑚, respectively, from left to right. The amount of injected charge increases as the applied voltage increased. Few amount of charges are observed inside the bulk with a field of -20 𝑘𝑉/𝑚𝑚, either because the net charge density is zero or because the detection techniques are restricted in sensitivity (0.1 𝐶/m 3 ). Positive and negative charges are detected immediately in the bulk after the voltage is applied at higher fields (i.e., -40, -60 𝑘𝑉/𝑚𝑚). These charges are injected at the anode and cathode, respectively, and move through the dielectric towards the opposite electrode.

Figure 37 :

 37 Figure 37: Space charge profiles as a function of the time and sample thickness, for fields of -20, -40 𝑎𝑛𝑑 60 𝑘𝑉/ 𝑚𝑚, from left to right, at room temperature. LDPE sample is coated by a gold layer on both sides. The color bar provides charge density scale in 𝐶/𝑚 3 .

Figure 38 :

 38 Figure 38: Net density of charge for -60 𝑘𝑉/𝑚𝑚 applied electric field at 2 periods of time, 100 s (blue) and 1000 s (red) Figure 39 represents the absolute current density with respect to time for applied fields -20, -40 𝑎𝑛𝑑 -60 𝑘𝑉/𝑚𝑚. It is obvious that the current density increases with respect to the applied electric field. Besides, even after three hours of polarization, the experimental current density seems to be decreasing with respect to polarization time and does not attain a stationary state. Furthermore, all the experimental current density curves have nearly the same slope.

Figure 39 :

 39 Figure 39: Absolute current density versus time for applied fields of (blue): -20 𝑘𝑉/𝑚𝑚, (red): -40 𝑘𝑉/𝑚𝑚, (yellow): -60 𝑘𝑉/𝑚𝑚 IV.4 Problem formulation IV.4.1 Selection of model parametersAll the unknown parameters of the BCT defined in Chapter 1 will be optimized in this section.

Figure 40 :

 40 Figure 40: Experimental protocol showing the initial net density of charges before each experiment.

  is ten decades lower than the experimental value. This theoretical value does not account the local interface conditions (local field strengthening on specific locations, chemical impurities forming deep traps) and therefore cannot be used in simulations. Other research proposed different values of the injection barriers in the range [1 𝑒𝑉 -1.

Figure 41 :

 41 Figure 41: Principle of the optimization technique. 𝑃 is the set of unknown parameters, 𝑥 represents the position, 𝑡 represents the time, 𝑓 𝑒𝑥𝑝 are the experimental data, 𝑓 𝑠𝑖𝑚 are the simulated data.

Figure 42 :

 42 Figure 42: Good and bad starting points of optimization

Figure 43 :

 43 Figure 43: The strategy used for finding the unique optimal solution that fit all experiments. 𝐸𝐶 𝑖=1..3 are three different experimental conditions, 𝑇 𝑖 is the temperature, 𝐸 𝑖 is the applied electric field, 𝑡 𝑝 𝑖 is polarization time

Figure 44 :

 44 Figure 44: Experimental current density (blue) versus simulated current density (red) using the optimal parameters 𝑃 1 , for applied fields of -20, -40, -60 𝑘𝑉/𝑚𝑚, from left to right.

  Figure 45 represents the space charge profiles obtained from simulating the three applied fields in (-20, -40, -60 𝑘𝑉/𝑚𝑚) using the BCT model using 𝑃 1 .

Figure 45 :

 45 Figure 45: Space charge profiles as a function of the time and sample thickness at room temperature, for fields of -20, -40 𝑎𝑛𝑑 -60 𝑘𝑉/𝑚𝑚, from left to right. Upper row: experimental charge density (PEA). Lower row: simulated charge density using 𝑃 1

Figure 46 :

 46 Figure 46: Space charge profiles as a function of the time and sample thickness at room temperature, for fields of -20, -40 𝑎𝑛𝑑 -60 𝑘𝑉/𝑚𝑚, from left to right. Upper row: experimental charge density (PEA). Lower row: simulated charge density using 𝑃 2

Figure 47 :

 47 Figure 47: Experimental current density (blue) versus simulated current density (red) using the optimal parameters 𝑃 2 , for applied fields of -20, -40, -60 𝑘𝑉/𝑚𝑚, from left to right.

Chapter IV 127 𝑃 3

 1273 are represented in Figure 48 and Figure 49. All the results are discussed in the following sections.

Figure 48 :

 48 Figure 48: Experimental current density (blue) versus simulated current density (red) using the optimal parameters 𝑃 3 , for applied fields of -20, -40, -60 𝑘𝑉/𝑚𝑚 (from right to left).

Figure 49 :

 49 Figure 49: Space charge profiles as a function of the time and sample thickness at room temperature, for fields of -20, -40 𝑎𝑛𝑑 -60 𝑘𝑉/𝑚𝑚, from left to right. Upper row: experimental charge density (PEA). Lower row: simulated charge density using 𝑃 3 It is clear from Figure 49 that 𝑃 3 can also produce a good match between experimental and simulated charge density for all of the three applied fields in (-20, -40, -60 𝑘𝑉/𝑚𝑚). Most of the properties highlighted in the experiment, such as the dominance of negative charge at a short time, then the dominance of positive charges until the end of the experiment, are reflected by the model.

Figure 50 :

 50 Figure 50: recombination rate between different types of carriers.

Figure 50

 50 Figure 50 represents the recombination rate between different kinds of carriers with respect to the applied electric field. The outputs are simulated using the BCT model using 𝑃 3 . The blue curve represents the recombination rate between trapped electrons and mobile holes; the red curve represents the recombination rate between mobile electrons and mobile holes; and finally, the black curve represents the recombination between mobile electrons and trapped holes. The curves demonstrate that recombination between trapped electrons and mobile holes (blue) has the lowest recombination rate compared to other recombination processes. Furthermore, the recombination rate between mobile electrons and mobile holes (red) is equivalent to that between mobile electrons

  It is worth noting that these results are based on only three different applied fields (-20, -40, and -60 𝑘𝑉/𝑚𝑚), which may be inaccurate. Additional experiments should be indroduced to the optimization process to improve the estimation of trapping and detrapping rates behavior at different fields, which could contribute to construct an empirical law able to characterize the trapping and detrapping rates in terms of the electric field.

Figure 51 :IV. 5 . 2 . 3 Mobility

 51523 Figure 51: Trapping and detrapping rates in 𝑠 -1 for electrons and holes as a function of applied field at 20°CIV.5.2.3 MobilityThe electrons and holes mobility could be estimated using Figure52by estimating the velocity of negative (blue) and positive (red) charges at short time (𝑡 = 10 𝑠). The velocity could be experimentally estimated by computing the slope equation which is found by determining the

Figure 52 :

 52 Figure 52: Experimental net density of charge. Applied field: -20 𝑘𝑉/𝑚𝑚. First 20 seconds of polarization process, temperature: 20 °C, sample thickness: 150 µ𝑚, color bar provides charge density scale in 𝐶/𝑚 3 . 𝑥 0 = 0 µ𝑚 for positive charges at 𝑡 = 0 and 𝑥 0 = 150 µ𝑚 for negative charges at 𝑡 = 0. 𝑥 1 represents the position at 𝑡 = 10 𝑠.

For

  the electrons: 𝑣 𝑒 (-20𝑘𝑉𝑚𝑚 -1 ) ≈ 7.2 𝜇𝑚/𝑠 For the holes: 𝑣 ℎ (-20𝑘𝑉𝑚𝑚 -1 ) ≈ 2.8 𝜇𝑚/𝑠 where 𝑣 𝑒 and 𝑣 ℎ are the approximated experimental velocities of electrons and holes respectively Chapter IV 132 Now to compute the experimental mobility of the electrons and holes, one has to substitute the velocity values in Equation 41: 𝑣 𝑒,ℎ = 𝜇 𝑒,ℎ 𝑒𝑥𝑝 × |𝐸 0 | (41)

  were approximated by calculating the slope of charge movement over a short time (𝑡 ≤ 10 𝑠), then these values were compared to the values provided by the BCT model. The experimental and simulated mobility values were found quite near to each other and varied in the same range. The outcomes of mobility values could be used as a starting point of optimization process in our future work.

Figure 53 .

 53 Figure 53. The experimental protocol aims to measure the charge density in LDPE of thickness 150 µm at two levels of fields, of -40 and -60 𝑘𝑉/𝑚𝑚, using the PEA method. As shown in Figure 53, each measurement consists of a polarization time of 90 minutes and a short-circuiting period of 24 hours between them.

Figure 53 :

 53 Figure 53: Material preparation and experimental protocol

Figure 54 :

 54 Figure 54: Space charge profiles as a function of the time and sample thickness, for fields of -40 𝑎𝑛𝑑 -60 𝑘𝑉/𝑚𝑚, from left to right, at room temperature. LDPE sample is coated by a gold layer on both sides.

Figure 55 :

 55 Figure 55: Space charge profiles as a function of the time and sample thickness, for fields of -40 and -60 𝑘𝑉/𝑚𝑚, from left to right, at room temperature. LDPE sample is coated by an Aluminum layer on both sides. According to Figure 55, for -40 𝑘𝑉/𝑚𝑚 (Left) applied field, negative charges are observed immediately after the voltage is applied. Only a few positive charges are injected compared to the quantity of negative charges through aluminum surfaces. Negative charges have the advantage inside the bulk of insulation compared to positive charges at any time of polarization. Moreover,

Figure 56 :

 56 Figure 56 : Space charge profiles as a function of the time and sample thickness, for fields of -40 and -60 𝑘𝑉/ 𝑚𝑚, from left to right, at room temperature. LDPE sample is coated by a copper layer on both sides.

Figure 57 :

 57 Figure 57: Space charge profiles as a function of the time and sample thickness, for fields of -40 and -60 𝑘𝑉/𝑚𝑚, from left to right, at room temperature. LDPE sample is coated by a Semi-conductor layer on both sides.

Figure 58 :

 58 Figure 58: Experimental protocol showing the initial net density of charges before each experiment

  min 𝑃∈ℝ 𝑛 𝐶(𝑃) 𝑠. 𝑡. 𝑙𝑏 ≤ 𝑃 ≤ 𝑢𝑏

Figure 59 :

 59 Figure 59: Space charge profiles as a function of the time and thickness of the sample, for fields of -40 𝑘𝑉/𝑚𝑚 (left column) and -60 𝑘𝑉/𝑚𝑚 (right column) at room temperature. Upper row: experimental measurements

Figure 60 :

 60 Figure 60: Space charge profiles as a function of the time and thickness of the sample, for fields of -40 𝑘𝑉/𝑚𝑚 (left column) and -60 𝑘𝑉/𝑚𝑚 (right column) at room temperature. Upper row: experimental measurements (PEA). Lower row: simulated results. The LDPE sample is coated by aluminum layer at both sides.

Figure 61 :

 61 Figure 61: Space charge profiles as a function of the time and thickness of the sample, for fields of -40 𝑘𝑉/𝑚𝑚 (left column) and -60 𝑘𝑉/𝑚𝑚 (right column) at room temperature. Upper row: experimental measurements

Figure 62 : 1

 621 Figure 62: Space charge profiles as a function of the time and thickness of the sample, for fields of -40 𝑘𝑉/𝑚𝑚 (left column) and -60 𝑘𝑉/𝑚𝑚 (right column) at room temperature. Upper row: experimental measurements (PEA). Lower row: simulated results. The LDPE sample is coated by semi-conductor layer at both sides. V.3.3 Discussion V.3.3.1 Space charge distribution According to Figure 59, 60, 61 and 62, 𝑃 ⋆ provided a good correlation between experimental and simulated charge density data for -40 𝑘𝑉/𝑚𝑚 and -60 𝑘𝑉/𝑚𝑚 for all considered surface

  concerned with understanding the transport mechanisms in solid dielectrics in order to predict their behavior under electrical stress. Experimental measurement alone does not clearly represent the microscopic phenomena inside dielectrics under electrical stress, such as trapping, detrapping, and recombination processes. Consequently, a BCT model was developed by our research team (DSF) that helps in predicting the dynamic behavior of LDPE under DC stress. This model provides a realistic representation of physical, electrical, and chemical phenomena at a microscopic scale that could not be accessible by experimental approaches. Our BCT model considers the trapping, de-trapping, injection, mobility, and recombination processes of positive and negative carriers. These carriers are generated by a modified Schottky injection law at the metal-dielectric interface and are extracted without any barriers. The BCT model is based on three basic equations that are known as the transport, Poisson, and continuity equations. Such models require some initial experimental conditions, such as temperature, applied voltage, polarization time, etc., alongside a set of parameters, such as injection barriers, mobility of carriers, trapping

  Experimental data that has been used includes net density charge measured by the Pulsed Electro-Acoustic method along with what is known as current measurements for measuring the current density. Moreover, the simulated data of charge and current densities were obtained by the BCT model developed for LDPE under DC stress. Five distinct optimization algorithms were discussed and compared before choosing the most suitable one to optimize the BCT model. The best algorithm was chosen by analyzing and comparing the robustness, accuracy, and efficiency of each algorithm, which reveals the strong Conclusion 158

Table 1 :

 1 Electrical and mechanical characteristics of low-density polyethylene[24]].

	Material	LDPE
	Density (g ⋅ cm -3 )	0.92
	Rate of crystallinity (%)	55-70
	Electrical rigidity at 23 °C (kV. mm -1 )	80
	Permittivity at 23 °C and 50Hz	2.3
	Glass transition temperature (°C)	-110
	Melting temperature (°C)	110-120

There are several types of polyethylene, which are differentiated by their chemical structure. Highdensity polyethylene (HDPE), obtained at low pressure, is a linear crystalline, rigid polymer where Chapter I 27 molecules can pack more closely. Low-density polyethylene (LDPE), produced at higher pressures, is a tough and flexible polymer characterized by long branches and semi-crystalline.

I.4.1 Basic equations I.4.1.1 Transport equation

	Charge transport is described by a hopping mechanism in which carriers move from site to site by
	tunneling a potential barrier. Mobile electrons and holes are associated with effective mobility.
	This mobility accounts for the possible trapping and de-trapping in shallow traps with a short-term
	residence. To introduce the hopping mechanism, the apparent effective mobilities 𝜇 𝑒,ℎ in equation
	(14) are defined by the depth of shallow traps 𝑤 𝑚𝑜 𝑒,ℎ :
	𝜇 𝑒,ℎ (𝑥, 𝑡) =	2𝑣𝑎 𝑒,ℎ 𝐸(𝑥, 𝑡)	• exp (-	𝑞𝑤 𝑚𝑜 𝑒,ℎ 𝑘 𝐵 𝑇	) • sinh ( 𝑞𝐸(𝑥, 𝑡)𝑎 𝑒,ℎ 2𝑘 𝐵 𝑇	)
									transport equation, poisson's
	equation (12), and continuity equation (13), below, considering a 1D problem along the spatial
	coordinate 𝑥.
	𝑗 𝑒,ℎ (𝑥, 𝑡) = 𝑞 ⋅ 𝜇 𝑒,ℎ (𝑥, 𝑡) • 𝐸(𝑥, 𝑡) • 𝑛 𝑒,ℎ (𝑥, 𝑡)	(11)
				1 𝑞		•	𝜕 2 𝑉(𝑥, 𝑡) 𝜕𝑥 2 =	-𝑛(𝑥, 𝑡) 𝜀 0 𝜀 𝑟	(12)
	𝜕𝑛 𝑒,ℎ (𝑥, 𝑡) 𝜕𝑡	+	1 𝑞	𝜕𝑗 𝑒,ℎ (𝑥, 𝑡) 𝜕𝑥	-	𝐷𝜕 2 𝑛 𝑒,ℎ (𝑥, 𝑡) 𝜕𝑥 2	= 𝑠 𝑖 (𝑥, 𝑡)	(13)

I.4.1.3 Continuity equations

									Chapter I
	𝑛(𝑥, 𝑡) = 𝑛 ℎ𝜇 (𝑥, 𝑡) + 𝑛 ℎ𝑡 (𝑥, 𝑡) -𝑛 𝑒𝜇 (𝑥, 𝑡) -𝑛 𝑒𝑡 (𝑥, 𝑡)	(18)
	Progressive change in charge density in time and space can be expressed as continuity equations
	for charge transport, such as:
	𝜕𝑛(𝑥, 𝑡) 𝜕𝑡	+	1 𝑞	𝜕𝑗(𝑥, 𝑡) 𝜕𝑥	-	𝐷𝜕 2 𝑛(𝑥, 𝑡) 𝜕𝑥 2	= 𝑠 𝑖 (𝑥, 𝑡)	(19)
	As mentioned previously, the diffusion phenomena is neglected in our model. Equation (19) will
	be substituted by Equation (20):
	𝜕𝑛 𝑒,ℎ (𝑥, 𝑡) 𝜕𝑡	+	1 𝑞	𝜕𝑗 𝑒,ℎ (𝑥, 𝑡) 𝜕𝑥	= 𝑠 𝑖 (𝑥, 𝑡)
							𝜕𝑉(𝑥,𝑡) 𝜕𝑥 . t is the time and x is the position.
	However, the net charge density is written as:
									38

  𝑛 ℎ𝜇 are the densities of mobile electrons and holes, 𝑛 𝑒𝑡 and 𝑛 ℎ𝑡 are the densities of trapped electrons and holes. Equation (21) can be applied to all species present in the material.

											Chapter I
	The de-trapping rate is defined by a de-trapping coefficient 𝑤 𝑡𝑟 𝑒,ℎ , such as:
									𝐷 𝑡𝑟 𝑒,ℎ = 𝑣 ⋅ exp (	-𝑤 𝑡𝑟 𝑒,ℎ 𝑘 𝐵 𝑇	)
		1 𝑞	⋅	𝑑𝑗 𝑒µ (𝑥, 𝑡) 𝜕𝑥	= 𝑠 𝑒𝜇 = -	𝑅 1 𝜇 𝑒 𝜀 0 𝜀 𝑟	𝑛 ℎ𝑡 𝑛 𝑒𝜇 -	𝑅 3 (𝜇 𝑒 + 𝜇 ℎ ) 𝜀 0 𝜀 𝑟	𝑛 𝑒𝜇 𝑛 ℎ𝜇 -𝐵 𝑒 𝑛 𝑒𝜇 (1 -	𝑛 𝑒𝑡 𝑁 0𝑒𝑡	) + 𝐷 𝑡𝑟 𝑒 𝑛 𝑒𝑡
	𝜕𝑛 𝑡𝑒 𝜕𝑡 𝜕𝑛 ℎµ 𝜕𝑡	+	1 𝑞	⋅	𝜕𝑗 ℎµ (𝑥, 𝑡) 𝜕𝑥	= 𝑠 𝑒𝑡 = -= 𝑠 ℎ𝜇 = -	𝑅 2 𝜇 ℎ 𝜀 0 𝜀 𝑟 𝑅 2 𝜇 ℎ 𝜀 0 𝜀 𝑟	𝑛 𝑒𝑡 𝑛 ℎ𝜇 -𝑅 0 𝑛 𝑒𝑡 𝑛 ℎ𝑡 + 𝐵 𝑒 𝑛 𝑒𝜇 (1 -𝑛 𝑒𝑡 𝑛 ℎ𝜇 -𝑅 3 (𝜇 𝑒 + 𝜇 ℎ ) 𝜀 0 𝜀 𝑟 𝑛 𝑒𝜇 𝑛 ℎ𝜇 -𝐵 ℎ 𝑛 ℎ𝜇 (1 -𝑛 𝑒𝑡 ) -𝐷 𝑡𝑟 𝑒 𝑛 𝑒𝑡 𝑁 0𝑒𝑡 𝑁 0ℎ𝑡 𝑛 ℎ𝑡	) + 𝐷 𝑡𝑟 ℎ 𝑛 ℎ𝑡	(21)
	𝜕𝑛 ℎ𝑡 𝜕𝑡						= 𝑠 ℎ𝑡 = -	𝑅 1 𝜇 𝑒 𝜀 0 𝜀 𝑟	𝑛 ℎ𝑡 𝑛 𝑒𝜇 -𝑅 0 𝑛 𝑒𝑡 𝑛 ℎ𝑡 + 𝐵 ℎ 𝑛 ℎ𝜇 (1 -	𝑛 ℎ𝑡 𝑁 0ℎ𝑡	) -𝐷 𝑡𝑟 ℎ 𝑛 ℎ𝑡
	𝑅 𝑖=0..3 are the reduction pre-factors multiplied by the Langevin's recombination rate, 𝐵 𝑒,ℎ are the
	trapping rates for electrons and holes, 𝐷 𝑡𝑟 𝑒,ℎ are the de-trapping rates for electrons and holes, 𝑁 0𝑒𝑡
	and 𝑁 0ℎ𝑡 are the maxima trapped charge densities for electrons and holes, respectively. 𝑛 𝑒𝜇 and
											39

Table 2 :

 2 Characteristics of the physical models for bipolar transport from 1994 to 2021. Model improvements by each researcher are written in bold.

	Reference	Charge Injection	Charge Extraction	Charge Transport	Charge Trapping	Charge recombination	Other	Parameters
	Alison [49]	Constant source at both electrodes	Non-blocking electrodes	Constant effective mobility	One deep trapping level, no detrapping	For mobile and trapped charges		Symmetric
	Fukuma [50]

Schottky injection at both electrodes Extraction barriers Hopping conduction between sites of the same energy

							Joule heating
					One deep trapping level, no detrapping	For mobile carriers	effects initial bulk accounted for	Symmetric
							charges
				Hopping
	Kaneko [51]	Schottky injection at both electrodes	Non-blocking electrodes	conduction between sites of the same	no deep trapping	For mobile carriers		Non-sym
				energy
	Le Roy [52]	Schottky injection at both electrodes	Non-blocking electrodes	Constant effective mobility	Trapping on detrapping one deep level with	For mobile and trapped charges	Initial bulk charges	Non-sym
				Hopping
	Boufayed [34]	Schottky injection at both electrodes	Non-blocking electrodes	mobility dependent on charge density	Exponential distribution of trap levels	For mobile and trapped charges		Non-sym
				and field
	Baudoin [53]	Schottky injection at both electrodes	Non-blocking electrodes	Constant effective mobility	Exponential distribution of trap levels

Langevin form for mobile and trapped charges Steady state Non-sym and sym

	Chen [64]	Schottky injection at both electrodes	Non-blocking electrodes	Constant effective mobility	Trapping on one deep level with detrapping	For mobile and trapped charges	corona-charging setup	Symmetric
				Hopping
	Baudoin [55] Injection rule	Non-blocking electrodes	mobility dependent on charge density	Exponential distribution of trap levels	For mobile and trapped charges	EL intensity	Symetrical
				and field
				Hopping
	Le Roy [35]	Modified Schottky at both electrodes	Non-blocking electrodes	mobility dependent on charge density	Exponential distribution of trap levels	For mobile and trapped charges
				and field

Temperature distribution equation (steady state) Non-Sym Zhan [56] Schottky at bottom and top + Threshold

	Alison and Hill (1994): In 1994, [49] Alison and Hill published a pioneering model to reproduce
	accumulations of space charges on a 2.5 mm sample of cross-linked polyethylene (XLPE), that
	considered transient procedures for charge generation and transport in cable insulation exposed to
	DC stresses (XLPE). The model included charging generation through injection into and through
	material bulk related to trapping and recombination at insulation-electrode interfaces.
	The principle for all model space charges lies in the explanation of the material charge conduction
	and mechanism of the electrical transport. The Alison-Hill model attempts to explain bipolar
	transport and spatial charging in solid dielectrics with high dc tension. The transport of bipolar
	materials is defined by an efficient mobility mechanism. This feature distinguishes the model from
	others. Charging carriers are effectively injected through electrodes, electrons from the cathode
	and holes from the anode. Injection takes the form of the Schottky mechanism, which overcomes
	a possible interface obstacle. After accessing the material, the carriers drift through the material
	characterized by effective mobility under the influence of the applied field. During the movements,
	many carriers are stuck in deep trap centers in localized states, thus reducing the overall number
	of charges. However, in the model there is no extraction barrier. On the other hand, they are
			Non-blocking electrodes	Hopping mobility dependent on charge density and field	Exponential distribution of trap levels	For mobile and trapped charges	heat transfer transient) equation (thermal	Non-Sym
	Doedens [57]	Modified Schottky and Fowler-Nordheim	Non-blocking electrodes	Hopping mobility dependent on charge density and field	Exponential depending on E distribution of trapping and detrapping levels. Trap rate	For mobile and trapped charges	Diffusion, trap three levels of	Non-Sym
	Zhan [58]	Schottky at bottom and top + Threshold	Non-blocking electrodes	Constant effective mobility	Exponential distribution of trap levels	For mobile and trapped charges	Ionization, Diffusion	Non-Sym

I 42 K. Kaneko et al. (1999): Kaneko

  [START_REF] Kaneko | Computer simulation on formation of space charge packets in XLPE films[END_REF] presented the macroscopic model with related assumptions to those of the model of Alison and Hill. The used material for simulation is LDPE. Bipolar charges

	are produced by a mechanism of the Schottky law and transported by hopping. There is no deep
	trap and the extraction barrier is assumed to be negligible. The recombination between holes and
	electrons is considered. The outcomes are space charge distribution and current density curve,
	compared to experimental measurements. The hypotheses did not allow them to achieve results
	that were consistent with the experimental data, emphasizing the importance of taking deep
	trapping into account.
	S

. Le Roy et al. (2004-2005): Le

	Roy [66] has proposed a model to characterize charge transport,
	trapping and recombination phenomena, experimentally observed by charge and discharge current
	measurements, space charge profiles and electroluminescence measurements, for LDPE. No
	extraction barrier has been applied.

. Boufayed et al. (2006): Instead

	interdependent processes (transportation, trapping, variation in mobility) and the complexities of
	the considered material
	F. Baudoin et al. (2008): F. Baudoin [53] has provided a method for a direct resolution in the
	steady-state of a bipolar charge transport model to understand how the various processes involved
	in the model act on the simulated current-voltage characteristic. Charge generation, transfer,
	trapping, and recombination are all included in the model. The steady-state solution obtained using
	of two single trap levels (shallow and low traps) used in other
	simulations, a model was developed by Boufayed et al. [63], [64], which introduced a more
	practical exponential distribution of traps. The study material was XLPE. The distribution
	parameters are based on distribution of Quirke et al. [65], [66] and molecular simulation of the
	physical and chemical traps. These traps are discreetly distributed and transformed into continuous
	distribution. A hopping mechanism is used for describing the transport. Quirke et al. assume that
	the traps are filled up from the lower energy level. Here, the recombination is not considered. The
	simulated results do not quite correlate with the experiment [34], which is explained by very

  Particular forms of the objective function and the constraints give rise to specialized algorithms that are more efficient. From this point of view, five different types

	of optimization models are discussed in this section, that are, Unconstrained optimization problems
	(2.5.1.1), Constrained optimization problems (2.5.1.2), Nonlinear programming problems
	(2.5.2.1), Nonlinear programming problems (2.5.2.2), and Least Squares optimization problems
	(2.5.3).

  Find a step length 𝛼 𝑘 , such that 𝑓(𝑥 𝑘 + 𝛼 𝑘 𝑥 𝑘 ) < 𝑓(𝑥 𝑘 ) 𝑑 𝑘 is usually calculated by a standard method, such as gradient descent 𝑑 𝑘 = -𝛻𝑓(𝑥 𝑘 )

			Chapter II
	or newton's method 𝑑 𝑘 = -	∇𝑓(𝑥 𝑘 ) 𝛻 2 𝑓(𝑥 𝑘 )
	4:	Update the design variables: 𝑥 𝑘+1 = 𝑃 𝑘 + 𝛼 𝑘 𝑑 𝑘
	5: 𝑘 ← 𝑘 + 1
	6: return 𝑥 𝑘

Algorithm 1: General for unconstrained problems Input condition: 𝑥 0 1: while stop criteria is not achieved, do 2: Compute a search direction 𝑑 𝑘 such that 𝛻𝑓(𝑥) 𝑑 𝑘 < 0 3: 62 directions

  Calculate 𝛼 𝑘 the solution of 𝐽 𝑟 (𝑃 𝑘 ) ⊤ 𝐽 𝑟 (𝑃 𝑘 )𝛼 𝑘 = -𝐽(𝑃 𝑘 )𝑟(𝑃 𝑘 )

		Chapter II
	3:	𝑃 𝑘+1 ← 𝑃 𝑘 + 𝛼 𝑘
	4: 𝑘 ← 𝑘 + 1
	5: return 𝑃 𝑘
		67

1: while stop criteria is not achieved, do 2:

Algorithm 1 :

 1 Trust Region method

					Chapter III
	Initial conditions: 𝑃 0 ∈ ℝ 𝑛 , 𝜂 1 , 𝜂 2 , 𝑟 0
	1: while stop criteria is not achieved, do
	2: Calculate 𝑠 𝑘 ⋆ by solving the subproblem (3) by using Dogleg or Steigaug-Toint
	methods
	3: Calculate 𝑙 =	𝑓(𝑥 𝑘 +𝑠 𝑘 ⋆ )-𝑓(𝑥 𝑘 ) 𝑚 𝑘 (𝑥 𝑘 +𝑠 𝑘 ⋆ )-𝑚 𝑘 (𝑥 𝑘 )
	4: if 𝑙 < 𝜂 1 then
	5:	𝑥 𝑘+1 ← 𝑥 𝑘
	6:	𝑟 𝑘 ←	1 2	‖𝑠 𝑘 ⋆ ‖
	7: if not
	8:	𝑥 𝑘+1 ← 𝑥 𝑘 + 𝑠 𝑘
	9:	if 𝑙 ≥ 𝜂 2 , then
	10:	𝑟 𝑘+1 ← 2𝑟 𝑘
	11:	if not
	12:	𝑟 𝑘+1 ← 𝑟 𝑘
	13: 𝑘 ← 𝑘 + 1
	14: return 𝑥 𝑘
					81

  If 𝐽 𝑟 (𝑃 𝑘 ) ⊤ 𝐽 𝑟 (𝑃 𝑘 ) is singular or ill-conditioned, LM introduces an approximation of the hessian matrix which makes it always positive definite, such that, 𝐻 𝑘 = 𝐽 𝑟 (𝑃 𝑘 ) 𝑇 𝐽 𝑟 (𝑃 𝑘 ) + 𝜇 𝑘 𝐼 𝑛 , with 𝜇 𝑘 (damping parameter) chosen such that 𝐻 𝑘 is positive definite and 𝐼 𝑛 is an identity matrix. Then we solve, as usual, the Newton Step equation: 𝛻 2 𝐶(𝑃)𝑠 𝑘 = -𝛻𝐶(𝑃) ⟺ (𝐽 𝑟 (𝑃 𝑘 ) ⊤ 𝐽 𝑟 (𝑃 𝑘 ) + 𝜇 𝑘 𝐼 𝑛 ) 𝑠 𝑘 = -𝐽 𝑟 (𝑃 𝑘 )𝑟(𝑃 𝑘 )

	This approach is known as	the Levenberg-Marquardt	(LM) algorithm,
	which is detailed in the Algorithm 2.

  Calculate 𝜇 𝑘 such that (𝐽 𝑟 (𝑃 𝑘 ) 𝑇 𝐽 𝑟 (𝑃 𝑘 ) + 𝜇 𝑘 𝐼 𝑛 ) ≻ 0 8: solve (𝐽 𝑟 (𝑃 𝑘 ) ⊤ 𝐽 𝑟 (𝑃 𝑘 ) + 𝜇 𝑘 𝐼 𝑛 )𝑠 𝑘 = -𝐽 𝑟 (𝑃 𝑘 )𝑟(𝑃 𝑘 )

	9:	𝑃 𝑘+1 ← 𝑃 𝑘 + 𝛼 𝑘 𝑠 𝑘
	10: 𝑘 ← 𝑘 + 1
	11: return 𝑃 𝑘

1: while stop criteria is not achieved, do 2: Calculate 𝛼 𝑘 3: Factorization of 𝐽 𝑟 (𝑃 𝑘 ) ⊤ 𝐽 𝑟 (𝑃 𝑘 ) 4: if 𝐽 𝑟 (𝑃 𝑘 ) ⊤ 𝐽 𝑟 (𝑃 𝑘 ) ≻ 0 then 5: solve 𝐽 𝑟 (𝑃 𝑘 ) ⊤ 𝐽 𝑟 (𝑃 𝑘 )𝑠 𝑘 = -𝐽 𝑟 (𝑃 𝑘 )𝑟(𝑃 𝑘 ) 6: if not 7:

1 :

 1 while Stop criteria is not achieved, do 2: evaluate f at each of the initial points such that: 𝐶(𝑃 0 ) ≤ 𝐶(𝑃 1 ) ≤ ⋯ ≤ 𝐶(𝑃 𝑛 ) 3: Calculate the center of gravity 𝑃 𝑔 of the points 𝑃 0 , … , 𝑃 𝑛 4: Construct a symmetric point 𝑃 𝑠 = 𝑃 𝑔 + (𝑃 𝑔 -𝑃 𝑛 ) of 𝑃 𝑛 with respect to the center of gravity 𝑃 𝑔 .

				Chapter III
	5:	if 𝐶(𝑃 𝑠 ) < 𝐶(𝑃 0 ) then
	6:	Enlarge the simplex in the direction of 𝑃 𝑠 :
	7:	𝑃 𝑎 ← 𝑃 𝑔 + 2(𝑃 𝑔 -𝑃 𝑛 )
	8:	if 𝐶(𝑃 𝑎 ) < 𝐶(𝑃 𝑠 ) then
	9:	𝑃 𝑛 ← 𝑃 𝑎
	10:	if not
	11:	𝑃 𝑛 ← 𝑃 𝑠
	12:	if 𝐶(𝑃 𝑛-1 ) < 𝐶(𝑃 𝑠 ) then
	13:	Reduce the simplex:
	14:	𝑃 𝑟 ← 𝑃 𝑛 +	1 2	(𝑃 𝑔 -𝑃 𝑛 )
	15:	if 𝐶(𝑃 𝑟 ) < 𝐶(𝑃 𝑠 ) then
	16:	𝑃 𝑛 ← 𝑃 𝑟
	17: if not
	18:	𝑃 𝑛 ← 𝑃 𝑠
	19: return 𝑃 𝑛
				87

Table 3 :

 3 Units, Symbols and bounds of the unknown parameters.

	Parameters	Barrier height of injection	Trapping depth	Trapping coefficient	Detrapping coefficient
	Unit	eV		eV		s -1		s -1
	Symbol	𝑤 𝑒	𝑤 ℎ	𝑤 𝑚𝑜 𝑒	𝑤 𝑚𝑜 ℎ	𝐵 𝑒	𝐵 ℎ	𝐷 𝑒	𝐷 ℎ
	Lower Bound (lb)	1		0.3	0		0
	Upper Bound (ub)	1.3		0.72	1		1

Table 4 :

 4 Optimized parameters using five optimization algorithms. 𝑃 0 is the starting point of the optimization algorithm, 𝑃 𝐸𝑥𝑎𝑐𝑡 represents the parameters we are trying to find. 𝑃 𝑇𝑅𝑅𝐴 , 𝑃 𝐿𝑀 , 𝑃 𝑁𝑀 , 𝑃 𝐺𝐴 , and 𝑃 𝑃𝑆 , are the optimal parameters produced respectively by, TRRA, LM, NM, GA, and PS.

	parameters	𝑤 𝑒	𝑤 ℎ	𝑤 𝑚𝑜 𝑒	𝑤 𝑚𝑜 ℎ	𝑤 𝑡𝑟 𝑒	𝑤 𝑡𝑟 ℎ	𝐵 𝑒	𝐵 ℎ
	𝑃 𝑒𝑥𝑎𝑐𝑡	1.21	1.10	0.71	0.60	0.97	0.76	0.006	0.0030
	𝑃 0	1.15	1.15	0.65	0.65	0.84	0.84	0.001	0.0010
	Algorithms
	𝑃 𝑇𝑅𝑅𝐴	1.21	1.10	0.71	0.60	0.97	0.76	0.006	0.0001
	𝑃 𝐿𝑀	1.21	1.10	0.71	0.60	0.97	0.98	0.006	0.0000
	𝑃 𝑁𝑀	1.19	1.17	0.66	0.64	0.85	0.88	0.999	0.1523
	𝑃 𝐺𝐴	1.23	1.10	0.68	0.60	0.80	0.80	0.981	0.0702
	𝑃 𝑃𝑆	1.20	1.09	0.66	0.60	0.93	0.76	0.078	0.7501

Table 5 :

 5 The time of computation, the number of iterations, and the number of evaluated functions required by each algorithm to find the optimal solution. 𝐶(𝑃 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 ) represents the value of cost function at the optimal solution found by each algorithm.

	Algorithms	Elapsed Time (s)	Number of	Func-count	Relative error
			iterations		𝐶(𝑃 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 )
	TRRA	1435	73	666	8.2322 × 10 -9
	LM	403	21	205	0.0005
	NM	1196	1028	1601	0.3336
	GA	135185	339	68000	0.0119
	PS	530248	1600	128080	43.6800

Table 6 :

 6 NM outcomes using 2 different initial starting points

	parameters	𝑤 𝑒	𝑤 ℎ	𝑤 𝑚𝑜 𝑒	𝑤 𝑚𝑜 ℎ	𝑤 𝑡𝑟 𝑒	𝑤 𝑡𝑟 ℎ	𝐵 𝑒	𝐵 ℎ
	𝑃 𝑒𝑥𝑎𝑐𝑡	1.21	1.10	0.71	0.60	0.97	0.76	0.006 0.003
	𝑃 0	1.19	1.13	0.68	0.62	0.95	0.77	0.005 0.002
	𝑃 𝑁𝑀	1.21	1.10	0.71	0.60	0.97	0.76	0.006 0.016
	𝑃 0	1.30	1.20	0.60	0.68	0.80	0.85	0.100 0.200
	𝑃 𝑁𝑀	1.29	1.28	0.60	0.70	0.81	0.88	0.053 0.173

Table 7 :

 7 Optimized parameters using five optimization algorithms.𝑃 0 is the starting point of the optimization algorithm, 𝑃 𝐸𝑥𝑎𝑐𝑡 represents the parameters we are trying to find. 𝑃 𝑇𝑅𝑅𝐴 , 𝑃 𝐿𝑀 , 𝑃 𝑁𝑀 , 𝑃 𝐺𝐴 , and 𝑃 𝑃𝑆 , are the optimal parameters produced respectively.

	Parameters	𝑤 𝑒	𝑤 ℎ	𝑤 𝑚𝑜 𝑒	𝑤 𝑚𝑜 ℎ	𝑤 𝑡𝑟 𝑒	𝑤 𝑡𝑟 ℎ	𝐵 𝑒	𝐵 ℎ
	𝑃 0	1.15	1.15	0.65	0.65	0.84	0.84	0.001	0.001
	Algorithms
	𝑃 𝑇𝑅𝑅𝐴	1.20	1.09	0.70	0.60	0.94	0.76	0.013	0.708
	𝑃 𝐿𝑀	1.19	1.10	0.59	0.61	1.17	0.66	0.207	-0.089
	𝑃 𝑁𝑀	1.22	1.13	0.71	0.6	0.84	0.89	0.154	0.080
	𝑃 𝐺𝐴	1.22	1.16	0.67	0.62	0.88	0.84	0.515	0.627
	𝑃 𝑃𝑆	1.20	1.09	0.60	0.61	0.93	0.77	0.999	0.031

Table 8 :

 8 The time of computation, the number of iterations, and the number of evaluated functions required by each algorithm to find the optimal solution. 𝐶(𝑃 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 ) represents the value of cost function at the optimal solution found by each algorithm.

	Algorithms	Computation	Number of	Func-count	Relative error
		Time (s)	iterations		𝐶(𝑃 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 )
	TRRA	1486	57	580	0.083
	LM	5820	26	289	0.353
	NM	3945	996	1554	0.111
	GA	154908	457	91600	0.133
	PS	688911	1761	158580	0.101

Table 9 :

 9 Comparison of algorithms with respect to robustness, efficiency, and accuracy

	Algorithm	Robustness	Efficiency	Accuracy
	TRRA

Table 10 :

 10 Units, Symbols and bounds of the unknown parameters.

	Parameters	Barrier height of injection	Shallow trapping depth	Trapping coefficient	Detrapping coefficient	Recombination reducation pre-factor	Inter-trap distance	Initial charge density
	Unit	eV	eV	s -1	s -1		unitless		𝑛𝑚	𝐶/𝑚 3
	Symbol	𝑤 𝑒 𝑤 ℎ	𝑤 𝑚𝑜 𝑒 𝑤 𝑚𝑜 ℎ	𝐵 𝑒 𝐵 ℎ	𝐷 𝑒 𝐷 ℎ	𝑅 1	𝑅 2	𝑅 3	𝑎 𝑒 𝑎 ℎ	𝛼
	Lower Bound (lb)	1	0.3	0	0		0.001		0	0
	Upper Bound (ub)	1.3	0.72	1	1		1		20	15

Table 11 :

 11 Trapping and detrapping parameters for different electric fields

	Electric field Trapping rate of	Trapping rate of	Detrapping rate of	Detrapping rate of
		electrons	holes	electrons	holes
	-20 𝑘𝑉/𝑚𝑚	𝐵 𝑒 -20	𝐵 ℎ -20	𝐷 𝑒 -20	𝐷 ℎ -20
	-40 𝑘𝑉/𝑚𝑚	𝐵 𝑒 -40	𝐵 ℎ -40	𝐷 𝑒 -40	𝐷 ℎ -40
	-60 𝑘𝑉/𝑚𝑚	𝐵 𝑒 -60	𝐵 ℎ -60	𝐷 𝑒 -60	𝐷 ℎ -60

Table 12 :

 12 Optimized parameters 𝑃 3 obtained by minimizing 𝐶 3 (𝑃)

	Parameters	𝑤 𝑒	𝑤 ℎ	𝑤 𝑚𝑜 𝑒	𝑤 𝑚𝑜 ℎ	𝑅 1	𝑅 2	𝑅 3	𝑎 𝑒	𝑎 ℎ
	𝑃 3	1.184	1.200	0.563	0.602	1	1	1	2.4	1.1

Table 13 :

 13 Optimized trapping and detrapping parameters obtained by minimizing 𝐶 3 (𝑃)

	Electric field	𝐵 𝑒	𝐵 ℎ	𝐷 𝑒	𝐷 ℎ
	-20 𝑘𝑉/𝑚𝑚	0.6770	0.0125	0.0030	0.0003
	-40 𝑘𝑉/𝑚𝑚	0.6772	0.0301	0.0034	0.0012
	-60 𝑘𝑉/𝑚𝑚	0.6802	0.0963	0.0035	0.0092

Table 14 :

 14 Optimal values of the initial net density of charges parameters by minimizing 𝐶 3 (𝑃)

	Parameters	𝛼 𝐸𝑥𝑝1 -60	𝛼 𝐸𝑥𝑝1 -40	𝛼 𝐸𝑥𝑝1 -20	𝛼 𝐸𝑥𝑝2 -60	𝛼 𝐸𝑥𝑝2 -40	𝛼 𝐸𝑥𝑝2 -20
	𝑃 3	0.01	0.02	0.03	1.25	5.49	14.7

  3 could be evaluated by substituting 𝑤 𝑚𝑜𝑒 , 𝑤 𝑚𝑜ℎ , 𝑎 𝑒 and 𝑎 ℎ in mobility equation, defined as:

	𝜇 𝑒,ℎ (𝑥, 𝑡) =	2𝑣𝒂 𝒆,𝒉 𝐸(𝑥, 𝑡)	• exp (-	𝑞𝒘 𝒎𝒐 𝒆,𝒉 𝑘 𝐵 𝑇	) • sinh ( 𝑞𝐸(𝑥, 𝑡)𝒂 𝒆,𝒉 2𝑘 𝐵 𝑇	)

Table 15 :

 15 Comparison between experimental and simulated mobility values of electrons and holes in 𝒎 𝟐 . 𝑽 -𝟏 . 𝒔 -𝟏 .

	Electric field	𝜇 𝑒	𝜇 𝑒 𝑒𝑥𝑝	𝜇 ℎ	𝜇 ℎ 𝑒𝑥𝑝
	-20 𝑘𝑉/𝑚𝑚	3.5 × 10 -13	3.6 × 10 -13	1.35 × 10 -14	1.4 × 10 -13
	-40 𝑘𝑉/𝑚𝑚	5 × 10 -13		1.5 × 10 -14
	-60 𝑘𝑉/𝑚𝑚	9 × 10 -13		1.8 × 10 -14

  ). It is worth noting that Experiment 1 started with highest applied electric field -60 𝑘𝑉/𝑚𝑚. The same for the current measurements (i.e., Experiment 2, Figure40), the initial charge density is found to be increasing after performing a

	new experiment (𝛼 𝐸𝑥𝑝2 -60 < 𝛼 𝐸𝑥𝑝2 -40 < 𝛼 𝐸𝑥𝑝2 -20 ). It should be mentioned that current measurements
	were performed after finishing the space charge measurements.
	The quantity of charges obtained during the current measurements are very high compared to the
	charges	found	initially	during	the	space	charge	measurements
	(i.e., 𝛼 𝐸𝑥𝑝2 -20 ≫ 𝛼 𝐸𝑥𝑝1

Table 16 ,

 16 17, 18 and 19.

	Chapter V

Table 16 :

 16 Units, Symbols, and range of the unknown parameters

	Parameters	Trapping depth	Recombination reduction pre-factors	Inter-Trap
	Unit	𝑒𝑉		Unitless		𝑛𝑚
	Symbol	𝑤 𝑚𝑜 𝑒 𝑤 𝑚𝑜 ℎ	𝑅 1	𝑅 2	𝑅 3	𝑎 𝑒 𝑎 ℎ
	Range	[0.3 -0.73]		[0.001 -1]	[0 -1]

Table 17 :

 17 Unknown injection barrier heights parameters for each surface electrode

	Electrode	Gold (Au)	Aluminum (Al)	Copper (Cu)	Semi-conductor (Sc)
	Symbols	𝑤 𝑒 𝐴𝑢 𝑤 ℎ 𝐴𝑢	𝑤 𝑒 𝐴𝑙	𝑤 ℎ 𝐴𝑙	𝑤 𝑒 𝐶𝑢	𝑤 ℎ 𝐶𝑢	𝑤 𝑒 𝑆𝑐	𝑤 ℎ 𝑆𝑐
	Unit	𝑒𝑉	𝑒𝑉		𝑒𝑉		𝑒𝑉
	Range	[0.9 -1.3]	[0.9 -1.3]	[0.9 -1.3]	[0.9 -1.3]

Table 18 :

 18 Unknown trapping and detrapping parameters for different electric fields

	Electric field	Trapping coefficient	Detrapping coefficient
	-40 𝑘𝑉/𝑚𝑚	𝐵 𝑒 -40 , 𝐵 ℎ -40	𝐷 𝑒 -40 , 𝐷 ℎ -40
	-60 𝑘𝑉/𝑚𝑚	𝐵 𝑒 -60 , 𝐵 ℎ -60	𝐷 𝑒 -60 , 𝐷 ℎ -60

Table 19 :

 19 Initial net charge density parameters for all types of surface electrodes

	Parameters	𝛼 𝐴𝑢	𝛼 𝐴𝑙	𝛼 𝐶𝑢	𝛼 𝑆𝑐
	Unit	𝐶/𝑚 3	𝐶/𝑚 3	𝐶/𝑚 3	𝐶/𝑚 3
	Range	[0 -15]	[0 -15]	[0 -15]	[0 -15]

Table 20 :

 20 Optimal values of the injection barrier height parameters for different electrodes

	Electrodes	Au	Al	Cu	Sc
	Parameters	𝑤 𝑒 𝐴𝑢 𝑤 ℎ 𝐴𝑢	𝑤 𝑒 𝐴𝑙 𝑤 ℎ 𝐴𝑙	𝑤 𝑒 𝐶𝑢 𝑤 ℎ 𝐶𝑢	𝑤 𝑒 𝑆𝑐 𝑤 ℎ 𝑆𝑐
	𝑃 ⋆	1.133 1.170	1.140 1.194	1.116 1.132	0.980 1.000

Table 21 :

 21 Optimal values of the initial net density of charges parameters

	Parameters	𝛼 𝐴𝑢	𝛼 𝐴𝑙	𝛼 𝐶𝑢	𝛼 𝑆𝑐
	𝑃 ⋆	0.1	0.2	0.1	4.7

Table 22 :

 22 Optimized parameters 𝑷 ⋆ obtained by minimizing 𝑪(𝑷)

	Parameters	𝑤 𝑚𝑜 𝑒 𝑤 𝑚𝑜 ℎ	𝑅 1	𝑅 2	𝑅 3	𝑎 𝑒	𝑎 ℎ
	𝑃 ⋆	0.71	0.62	0.001	1	1	5.08 1.88

Table 23 :

 23 Optimized trapping and detrapping parameters using current density measurements data.

	Electric field	𝐵 𝑒	𝐵 ℎ	𝐷 𝑒	𝐷 ℎ
	-40 𝑘𝑉/𝑚𝑚	0.0006	0.0279	0.0004	0.0026

If 𝐶(𝑃 1 ) ≤ 𝐶(𝑃 𝑠 ) < 𝐶(𝑃 2 ), reduce the simplex (outside), such that 𝑃 𝑟 = 𝑃 𝑔 + 1 2 (𝑃 𝑠 -𝑃 𝑔 ), as shown in Figure 20.