2. Thermoelectric effects

voltage is applied across a thermoelectric material. Seebeck and Peltier effects depend on each other, and this dependence was demonstrated by W. Thomson who also showed the existence of a third thermoelectric effect, known as the Thomson effect. Thomson effect describes reversible heating or cooling, in a homogeneous semiconductor material, when there is both a flow of electric current and a temperature gradient [2, 3]. For thermoelectric cooling devices (TECs), a thermocouple consists of a p-type and n-type legs, with Seebeck coefficients (α) values positives and negatives respectively, joined by a conductor metal with low α value; in this chapter, we take this value as zero for calculations. Practical devices make use of modules that contain many thermocouples connected electrically in series and thermally in parallel [4]. TECs suffer from low efficiency, therefore, research on system geometry, for design and fabrication of thermoelectric cooling devices, is investigated in recent days [5, 6]. Coefficient of performance (COP) is the most important parameter for a thermoelectric cooling device, which is defined as the heat extracted from the source due an electrical energy applied [7]. Single-stage devices operate between a heat source and sink at a temperature gradient. However, multistage devices provide an alternative for extending the maximum temperature difference for a thermoelectric cooler. Therefore, two-stage coolers should be used to improve the cooling power, Qc, and COP. In recent days, multistage thermoelectric cooling devices have been developed as many as six stages with bismuth telluride-based alloys. Recent works have investigated the ratio of the TE couple number between the stages and the effects of thermocouple physical size and have found that the cooling capacity is closely related to its geometric structure and operating conditions [8, 9]. In this chapter, a thermodynamics analysis and optimization procedure on performance of two-stage thermoelectric cooling devices based on the properties of established materials, system geometry and energy conversion, is analysed. Energy conversion issues in thermoelectric devices can be solved according to material properties: by increasing the magnitude of the differential Seebeck coefficient, by increasing the electrical conductivities of the two branches, and by reducing their thermal conductivities [10]. Several new theoretical and practical methods for the improvement of materials have been put forward and, at last, it seems that significant advances are being made, at least on a laboratory scale. In this work, we consider temperature-dependent properties material (TDPM) systems in calculations to determine the influence of the Thomson effect on performance [11, 12]. Many investigations have been conducted to improve the cooling capacity of two-stage TEC and found that cooling capacity is closely related to geometric structure and operating conditions of TECs. Our analysis to optimize cooling power of a thermoelectric micro-cooler (TEMC) includes a geometric optimization, that is, different cross-sectional areas for the p-type and ntype legs in both stages [13]. We find a novel procedure based on optimal material configurations, using two different semiconductors with different material properties, to improve the

290 Bringing Thermoelectricity into Reality

performance of a TEMC device with low-cost production.

This chapter is organized as follows: in Section 2, we give an overview of the thermoelectric effects. In Section 3, we apply thermodynamics theory to solve thermoelectric systems, and consequently, a description of the operation of a TEC device is presented. In Section 4, we proposed a two-stage TEC model taken into account Thomson effect for calculations to show its impact on COP and Qc. In Section 5, geometric parameters, cross-sectional area (A), and length (L) of a proposed two-stage TEMC system is analysed. For this purpose, constant Thermoelectricity results from the coupling of Ohm's law and Fourier's law. Thermoelectric effects in a system occur as the result of the mutual interference of two irreversible processes occurring simultaneously in the system, namely heat transport and charge carrier transport [14]. To define Seebeck and Peltier coefficients, we refer to the basic thermocouple shown in Figure 1, which consists of a closed circuit of two different semiconductors. For a thermocouple composed of two different materials a and b, the voltage is given by:

$$V\_{ab} = \int\_{1}^{2} (\alpha\_b - \alpha\_a) dT \tag{1}$$

where the parameters αa and αb are the Seebeck coefficients for semiconductor materials a and b.

Figure 1. Single thermocouple model for a TEC system.

The differential Seebeck coefficient, under open-circuit conditions, is defined as the ratio of the voltage, V, to the temperature gradient,ΔT

$$
\alpha\_{ab} = \frac{V}{\Delta T} \tag{2}
$$

3. Thermoelectric refrigeration in nonequilibrium thermodynamics

charge transport relations, consistent with the Onsager theory [17], are

el is the electric current density, j

∇ ! � κ ∇ ! T þ j

κð Þ T ∂2 T ∂x2 þ

∇ � j

j

j

Theory of thermoelectric cooling is analysed according to out-of-equilibrium thermodynamics. Under isotropic conditions, when an electrical current density flows through the semiconductor material with a temperature gradient and steady-state condition, the heat transport and

Thermoelectric Cooling: The Thomson Effect in Hybrid Two-Stage Thermoelectric Cooler Systems with Different…

where, α is the Seebeck coefficient, T is the temperature, κ is the thermal conductivity, E is the

Equation (9) is the essential equation for thermoelectric phenomena. The governing equations

q ¼ j

el ¼ 0 and ∇ � j

2 r � T dα dT J ! � ∇ !

where r is the electrical resistivity ð Þ r ¼ 1=σ and J is the electric current density. In Equation (11), the first term describes the thermal conduction due to the temperature gradient. According to Fourier's law, the second term is the joule heating and the third term is the Thomson heat, both depending on the electric current density [18]. Now, from Equation (11),

Thermoelectric coolers make use of the Peltier effect which origin resides in the transport of heat by an electric current. For this analysis, we assume that thermal conductivity, electrical resistivity, and Seebeck coefficient are all independent of temperature, that is, CPM model [19], and the metal that connects the p-type with the n-type leg has a low α value, therefore it is considered as zero. We assume that there is zero thermal resistance between the ends of the branches and the heat source and sink. Thus, only electrical resistance is considered for the thermocouple legs, thereby, the thermocouple legs are the only paths to transfer heat between the source and sink, conduction via the ambient, convection, and radiation are ignored. These

� jT dα dT ∂T ∂x ¼ � j

For one-dimensional model, from Equations (8) and (9), we get for the heat flux

the equation that governs the system for one-dimensional steady state is given by:

∂T ∂x 2

dκ dT

el ¼ σE � σα∇T (8)

q ¼ αTjel � κ∇T (9)

2

q is the heat flux and σ is the electric conductivity.

el � E (10)

http://dx.doi.org/10.5772/intechopen.75440

293

T ¼ 0 (11)

σð Þ T (12)

framework

and

are

electric field, j

3.1. Cooling power

Electrons move through the n-type element towards the positive pole, attraction effect, while the negative pole of the voltage source repels them. Likewise, in the p-type semiconductor, the holes move to the negative potential of the voltage source, while positive potential acts as repel of the holes and they move in the contrary direction to the flow of electrons. As a result, in ptype semiconductors, α is positive and in n-type semiconductors, α is negative [15]. Peltier coefficient is equal to the rate of heating or cooling, Q, ratio at each junction to the electric current, I. The rate of heat exchange at the junction is

$$Q = \pi\_{ab}I \tag{3}$$

Peltier coefficient is regarded as positive if the junction at which the current enters is heated and the junction at which it leaves is cooled. When there is both an electric current and a temperature gradient, the gradient of heat flux in the system is given by

$$\frac{dQ}{d\mathbf{x}} = \pi I \frac{dT}{d\mathbf{x}}\tag{4}$$

where x is a spatial coordinate and T the temperature. Thomson coefficient, known as the effect of liberate or absorb heat due to an electric current flux through a semiconductor material in which exist a temperature gradient, is given by the Kelvin relation as follows

$$
\tau\_a - \tau\_b = T \frac{d\alpha\_{ab}}{dT} \tag{5}
$$

When Seebeck coefficient is considered independent of temperature, Thomson coefficient will not be taken into account in calculations, τ is zero.

#### 2.1. Thomson relations

Seebeck effect is a combination of the Peltier and Thomson effects [16]. The relationship between temperature, Peltier, and Seebeck coefficient is given by the next Thomson relation

$$
\pi\_{ab} = \alpha\_{ab} T \tag{6}
$$

These last effects have a relation to the Thomson coefficient, τ, given by

$$
\pi = T \frac{d\alpha}{dT} \tag{7}
$$

To develop an irreversible thermodynamics theory, Thomson's theory of thermoelectricity plays a remarkable role, because this theorem is the first attempt to develop such theory.


**2. Anatomy**

## **2.1 Lacrimal gland and excretory system**

#### **2.1.1 Lacrimal gland**

The main lacrimal gland is located in a shallow depression along the superior lateral orbit. There is fibroadipose tissue between the gland and the orbit. The gland is divided into 2 parts by a lateral expansion of the levator apeunorosis. An isthmus of glandular tissue occasionally exists between the palpebral lobe and the main orbital gland1.

Many accessory lacrimal glands can be found along the inner surface of the eyelids. A variable number of thin-walled excretory ducts, blood vessels, lymphatics, and nerves pass from the main orbital gland into these accessory lacrimal glands. The ducts continue downward, and about 12 of them empty into the conjunctival fornix approximately 5 mm above the superior margin of the upper tarsus. Because the lacrimal excretory ducts pass through the palpebral portion of the gland, biopsy of the lacrimal gland is usually performed on the main part to avoid sacrificing the ducts1.

The lacrimal glands are exocrine glands, and they produce a serous secretion. The body of each gland contains 2 cell types:

acinar cells, which line the lumen of the gland

Endoscopic Dacryocystorhinostomy 51

Fig. 2. Position of lacrimal sac as seen during endonasal visualization

Occasionally, probing may be necessary to achieve patency.

resembles adenoid tissue and has a rich venous plexus and elastic fibers.

**3. Etiologies and predisposing factors of lacrimal obstruction**

lid. In 90% of people, the two canaliculi join to make a common canaliculus before entering the lacrimal sac. The lacrimal sac is placed in an oval-shaped fossa measuring 15 mm in height and 10 mm in width. This fossa is bounded by anterior and posterior lacrimal crests which fuse at a suture line that crosses the lacrimal fossa in a vertical manner. The lacrimal sac opens into nasolacrimal canal which is formed by the maxillary, lacrimal, and inferior turbinate bones. The nasolacrimal duct passes through this osseous canal for approximately 12 mm. Then it turns into a membranous duct for 5 mm before entering the inferior meatus2. The duct orifice is often covered by Hasner's valve to prevent reflux of secretions. In about 30% of full-term neonates, the outlet of the nasolacrimal duct is closed for up to 6 months.

The lacrimal puncta and the canaliculi are lined with stratified squamous nonkeratinized epithelium that merges with the epithelium of the eyelid margins. Near the lacrimal sac, the epithelium changes into 2 layers: a superficial columnar layer and a deep, flattened cell layer. Goblet cells and occasional cilia are present. In the canaliculi, the substantia propria consists of collagenous connective tissue and elastic fibers. The wall of the lacrimal sac

Patients with obstruction of lacrimal system usually complain of excessive tearing or epiphora. When dacryocystitis occurs, purulent drainage or inflammation can be noticed in

 myoepithelial cells, which surround the parenchyma and are covered by a basement membrane

The lacrimal artery, a branch of the ophthalmic artery, supplies the gland. The lacrimal gland receives secretomotor cholinergic, vasoactive intestinal polypeptide (VIP)-ergic,and sympathetic nerve fibers in addition to a sensory innervation via the lacrimal nerve (CN V1). Cyclic adenosine monophosphate is the second messenger for VIP and β-adrenergic stimulation of the gland; cholinergic stimulation acts through an inositol 1,4,5-triphosphateactivated protein kinase C. The gland also contains α1-adrenergic receptors. Extremely complex, the gland's neuroanatomy governs both reflex and psychogenic stimulation1.

#### **2.1.2 Accessory glands**

The accessory lacrimal glands of Krause and Wolfring are located at the proximal lid borders or in the fornices and are cytologically identical to the main lacrimal gland, receiving a similar innervation. They account for about 10% of the total lacrimal secretory mass1.

#### **2.1.3 Lacrimal excretory system**

The lacrimal excretory (drainage) system includes the upper and lower puncta, the lacrimal canaliculi, the lacrimal sac, and the nasolacrimal duct. It is important to note that the first 2 mm of canaliculi are perpendicular to the lid margin but the distal 8 mm are parallel to the

Fig. 1. Anatomy of the left lacrimal apparatus

1.superior punctum 2.superior canaliculus 3.inferior punctum 4.inferior canaliculus 5.medial canthal ligament 6.common canaliculus 7.lacrimal sac 8. Lacrimal duct 9.middle turbinate 10.lacrimal bone 11.inferior turbinate 12.Hasner's valve

myoepithelial cells, which surround the parenchyma and are covered by a basement

The lacrimal artery, a branch of the ophthalmic artery, supplies the gland. The lacrimal gland receives secretomotor cholinergic, vasoactive intestinal polypeptide (VIP)-ergic,and sympathetic nerve fibers in addition to a sensory innervation via the lacrimal nerve (CN V1). Cyclic adenosine monophosphate is the second messenger for VIP and β-adrenergic stimulation of the gland; cholinergic stimulation acts through an inositol 1,4,5-triphosphateactivated protein kinase C. The gland also contains α1-adrenergic receptors. Extremely complex, the gland's neuroanatomy governs both reflex and psychogenic stimulation1.

The accessory lacrimal glands of Krause and Wolfring are located at the proximal lid borders or in the fornices and are cytologically identical to the main lacrimal gland, receiving a similar innervation. They account for about 10% of the total lacrimal secretory

The lacrimal excretory (drainage) system includes the upper and lower puncta, the lacrimal canaliculi, the lacrimal sac, and the nasolacrimal duct. It is important to note that the first 2 mm of canaliculi are perpendicular to the lid margin but the distal 8 mm are parallel to the

1.superior punctum 2.superior canaliculus 3.inferior punctum 4.inferior canaliculus 5.medial canthal ligament 6.common canaliculus 7.lacrimal sac 8. Lacrimal duct 9.middle

turbinate 10.lacrimal bone 11.inferior turbinate 12.Hasner's valve

membrane

**2.1.2 Accessory glands**

**2.1.3 Lacrimal excretory system**

Fig. 1. Anatomy of the left lacrimal apparatus

mass1.

Fig. 2. Position of lacrimal sac as seen during endonasal visualization

lid. In 90% of people, the two canaliculi join to make a common canaliculus before entering the lacrimal sac. The lacrimal sac is placed in an oval-shaped fossa measuring 15 mm in height and 10 mm in width. This fossa is bounded by anterior and posterior lacrimal crests which fuse at a suture line that crosses the lacrimal fossa in a vertical manner. The lacrimal sac opens into nasolacrimal canal which is formed by the maxillary, lacrimal, and inferior turbinate bones. The nasolacrimal duct passes through this osseous canal for approximately 12 mm. Then it turns into a membranous duct for 5 mm before entering the inferior meatus2. The duct orifice is often covered by Hasner's valve to prevent reflux of secretions. In about 30% of full-term neonates, the outlet of the nasolacrimal duct is closed for up to 6 months. Occasionally, probing may be necessary to achieve patency.

The lacrimal puncta and the canaliculi are lined with stratified squamous nonkeratinized epithelium that merges with the epithelium of the eyelid margins. Near the lacrimal sac, the epithelium changes into 2 layers: a superficial columnar layer and a deep, flattened cell layer. Goblet cells and occasional cilia are present. In the canaliculi, the substantia propria consists of collagenous connective tissue and elastic fibers. The wall of the lacrimal sac resembles adenoid tissue and has a rich venous plexus and elastic fibers.

#### **3. Etiologies and predisposing factors of lacrimal obstruction**

Patients with obstruction of lacrimal system usually complain of excessive tearing or epiphora. When dacryocystitis occurs, purulent drainage or inflammation can be noticed in

Endoscopic Dacryocystorhinostomy 53

Dacryoliths (lacrimal stones) or "calculi" of the nasolacrimal ducts were described by Cesoni in as early as 16704, and have been reported to occur in between 6 and 18% of patients with nasolacrimal duct obstruction who undergo dacryocystorhinostomy (DCR). Dacryoliths may occur in any part of the nasolacrimal system, albeit most commonly in the lacrimal sac. Several predisposing factors have been suggested, such as increased occurrence in females, patient age below 50 years, association with cigarette smoking and facial-sinonasal trauma, and increased frequency subsequent to previous occurrence of dacryocystitis. However, other studies have indicated increased frequency in males and patients aged above 50 years. Therefore, it seems that both genders are involved to nearly the same extent. Dacryoliths usually become symptomatic when they obstruct the nasolacrimal system. This can result in epiphora, acute dacryocystitis, protrusion of the lacrimal canthal region, and partial closure of the lacrimal passage (recognized during syringing by the ophthalmologist). Interestingly, dacryoliths occur more often in patients with partial and incomplete closure of the lacrimal passage (i.e., patients with epiphora despite patent lacrimal passages on syringing). Scanning electron microscopy has shown that dacryoliths are composed of lobes and lobules built on an amorphous core material5. Atomic absorption spectrophotometric investigations demonstrate that dacryoliths consist almost entirely of organic proteins and, to a much lesser extent, of inorganic material5. According to Lew et al., lacrimal fluid from patients with dacryoliths contains a reduced amount of lysozyme and a lower calcium concentration than normal lacrimal fluid. It is important to recognize that daryoliths are not calcified or composed of any other "hard" substances. Some stones reveal hyphae-like structures,

A comprehensive ophthalmologic examination is mandatory in the primary evaluation of every patient with lacrimal system obstruction. An examination with the slit lamp can reveal the normal or abnormal tear film over the conjunctiva and if the thickness of the tearfilm is more than usual, it can be a sign of lacrimal drainage system obstruction. In addition, the ocular surface, eyelid structures, visual acuity, extraocular motility, and visual

Gentle pressure over the sac produces reflux of mucopurulent material suggestive of lower

Irrigation test is another useful test in assessing patients. In this test, an appropriate lacrimal syringe is passed through the inferior lacrimal punctum and 2-5 ml of sterile distilled water is injected and pushed though the inferior canaliculus. If the water passes easily into the nose and the patient senses that, the patency of the system is confirmed. Otherwise, it is one of the most reliable signs of lacrimal system obstruction. Some authors recommend that after either external or endoscopic DCR, this test can be performed indicating the patency of the system. Nasal examination, especially nasal endoscopy, should be obligatory for every lacrimal obstruction patient. The examination of the lacrimal area with the nasal speculum and headlight provides only a poor view of this region and is not sufficient; Endoscopy provides a clear diagnostic look for nasal polyps, imporant anatomic variations, tumors, and other

although no fungi were recovered by culturing5.

field should be tested and documented before surgery.

pathological endonasal conditions such as septal deviation.

**4. Assessment of the patient**

sac obstruction (regurgitation test).

**4.1 Physical examination**

the medial canthal region. It is important to ask patients about any nasal airway obstruction, drainage or epistaxis, which may suggest intranasal causes of lacrimal obstruction, such as polyps or neoplasms. Sometimes, the nasolacrimal duct is injured secondary to prior sinus surgery, particularly a large maxillary antrostomy.

Lacrimal excretory system foreign bodies are rare but they can impair draining function and might be presented as epiphora, recurrent attacks of acute dacryocystitis and in some patients, chronic dacryocystitis3. Exogenous foreign bodies in most patients lodge in lacrimal sac or nasolacrimal duct after external manipulation3. Foreign bodies in some patients have endogenous origin and in the form of dacryoliths may lead to lacrimal flow obstruction3. In both forms, surgical removal of foreign bodies is necessary3. Classic surgical approach is external dacryocystorhinostomy (DCR) but in recent years with rapid improvement of endoscopic techniques intranasal approaches introduce themselves as an effective substitute for external DCR3. These approaches are helpful for preoperative diagnosis and effective for surgical removal of lacrimal foreign bodies.

We had an experience with foreign bodies which was published in Iranian journal of ophthalmology. In our case, the lacrimal sac foreign body was a piece of silicon tube that was used as a stent in previous external DCR. On retrospective enquiry we found that at the time of silicon tube removal, it was pulled forcefully out through the lacrimal canaliculi and when it was impacted at the given site it was cut and the remaining part could not be found through the nose. The presenting signs and symptoms of this case were completely similar to the failed DCR procedure, so she was referred to our department for more evaluation. Anterior rhinoscopy is normal in many of such cases, so we should emphasize on the role of the nasal endoscopy as a safe and rapid diagnostic method. In the nasal endoscopy, the condition of rhinostomy site can be evaluated and any foreign body, granulation tissue, scar formation or synechia between middle turbinate and lateral nasal wall can be found and appropriate treatment plan can be established3.

Fig. 3. Rhinostomy site with foreign body T: Middle turbinate, R: Previous rhinostomy site, F: Foreign body

the medial canthal region. It is important to ask patients about any nasal airway obstruction, drainage or epistaxis, which may suggest intranasal causes of lacrimal obstruction, such as polyps or neoplasms. Sometimes, the nasolacrimal duct is injured secondary to prior sinus

Lacrimal excretory system foreign bodies are rare but they can impair draining function and might be presented as epiphora, recurrent attacks of acute dacryocystitis and in some patients, chronic dacryocystitis3. Exogenous foreign bodies in most patients lodge in lacrimal sac or nasolacrimal duct after external manipulation3. Foreign bodies in some patients have endogenous origin and in the form of dacryoliths may lead to lacrimal flow obstruction3. In both forms, surgical removal of foreign bodies is necessary3. Classic surgical approach is external dacryocystorhinostomy (DCR) but in recent years with rapid improvement of endoscopic techniques intranasal approaches introduce themselves as an effective substitute for external DCR3. These approaches are helpful for preoperative

We had an experience with foreign bodies which was published in Iranian journal of ophthalmology. In our case, the lacrimal sac foreign body was a piece of silicon tube that was used as a stent in previous external DCR. On retrospective enquiry we found that at the time of silicon tube removal, it was pulled forcefully out through the lacrimal canaliculi and when it was impacted at the given site it was cut and the remaining part could not be found through the nose. The presenting signs and symptoms of this case were completely similar to the failed DCR procedure, so she was referred to our department for more evaluation. Anterior rhinoscopy is normal in many of such cases, so we should emphasize on the role of the nasal endoscopy as a safe and rapid diagnostic method. In the nasal endoscopy, the condition of rhinostomy site can be evaluated and any foreign body, granulation tissue, scar formation or synechia between middle turbinate and lateral nasal wall can be found and

Fig. 3. Rhinostomy site with foreign body T: Middle turbinate, R: Previous rhinostomy site,

surgery, particularly a large maxillary antrostomy.

appropriate treatment plan can be established3.

F: Foreign body

diagnosis and effective for surgical removal of lacrimal foreign bodies.

Dacryoliths (lacrimal stones) or "calculi" of the nasolacrimal ducts were described by Cesoni in as early as 16704, and have been reported to occur in between 6 and 18% of patients with nasolacrimal duct obstruction who undergo dacryocystorhinostomy (DCR). Dacryoliths may occur in any part of the nasolacrimal system, albeit most commonly in the lacrimal sac. Several predisposing factors have been suggested, such as increased occurrence in females, patient age below 50 years, association with cigarette smoking and facial-sinonasal trauma, and increased frequency subsequent to previous occurrence of dacryocystitis. However, other studies have indicated increased frequency in males and patients aged above 50 years. Therefore, it seems that both genders are involved to nearly the same extent. Dacryoliths usually become symptomatic when they obstruct the nasolacrimal system. This can result in epiphora, acute dacryocystitis, protrusion of the lacrimal canthal region, and partial closure of the lacrimal passage (recognized during syringing by the ophthalmologist). Interestingly, dacryoliths occur more often in patients with partial and incomplete closure of the lacrimal passage (i.e., patients with epiphora despite patent lacrimal passages on syringing).

Scanning electron microscopy has shown that dacryoliths are composed of lobes and lobules built on an amorphous core material5. Atomic absorption spectrophotometric investigations demonstrate that dacryoliths consist almost entirely of organic proteins and, to a much lesser extent, of inorganic material5. According to Lew et al., lacrimal fluid from patients with dacryoliths contains a reduced amount of lysozyme and a lower calcium concentration than normal lacrimal fluid. It is important to recognize that daryoliths are not calcified or composed of any other "hard" substances. Some stones reveal hyphae-like structures, although no fungi were recovered by culturing5.


*Article* **Controllability Evaluation of EV Charging Infrastructure Transformed from Gas Stations in Distribution Networks with Renewables**

#### **Shuang Gao 1,\*, Jianzhong Wu 2,\* and Bin Xu 3**

Received: 13 March 2019; Accepted: 22 April 2019; Published: 25 April 2019

**Abstract:** A considerable market share of electric vehicles (EVs) is expected in the near future, which leads to a transformation from gas stations to EV charging infrastructure for automobiles. EV charging stations will be integrated with the power grid to replace the fuel consumption at the gas stations for the same mobile needs. In order to evaluate the impact on distribution networks and the controllability of the charging load, the temporal and spatial distribution of the charging power is calculated by establishing mapping the relation between gas stations and charging facilities. Firstly, the arrival and parking period is quantified by applying queuing theory and defining membership function between EVs to parking lots. Secondly, the operational model of charging stations connected to the power distribution network is formulated, and the control variables and their boundaries are identified. Thirdly, an optimal control algorithm is proposed, which combines the configuration of charging stations and charging power regulation during the parking period of each individual EV. A two-stage hybrid optimization algorithm is developed to solve the reliability constrained optimal dispatch problem for EVs, with an EV aggregator installed at each charging station. Simulation results validate the proposed method in evaluating the controllability of EV charging infrastructure and the synergy effects between EV and renewable integration.

**Keywords:** vehicle-to-grid; EV charging infrastructure; optimal dispatch; oil-to-electricity transformation for automobiles

#### **1. Introduction**

The high penetration of electric vehicles and renewable energy sources in the electric power networks provides a promising solution to the energy crisis and environmental stress. If properly controlled, EVs can act as mobile energy storage to facilitate the integration of intermittent renewable power generation, further alleviating the dependency on fossil fuels [1]. The feasibility of an electrified transportation system for zero or low carbon emission in the near future has been studied [2,3]. China Energy Administration announced the goal to have 120 thousand EV charging stations and 4.8 million charging piles in the power network by 2020. With increasing market share of electric vehicles (EVs), the gas stations are expected to be replaced with EV charging facilities supplied by electric power systems. However, as the EV charging introduces additional load demand, new challenges will be placed on the power system, consisting of reliability and power quality issues [4]. The charging load distribution transformed from the fuel consumption at the gas stations needs to be precisely quantified to evaluate the impact on the power system. Then the possible solutions can be suggested, such as upgrading the existing power system or applying appropriate control techniques in order to maintain the reliable and economic operation of the power system.

As the fuel consumption at the gas stations is transformed to electric energy consumption at various charging facilities, the temporal and spatial distribution of energy demand in the electric power network needs to be quantified. The economic benefits have been studied by using a lifetime cycle model for oil-to-electricity transformation in the transportation sector across the city [5,6]. However, the mapping relation between gas station and charging station for computing the distribution of charging load in the power network, based on the accessible service data of the facilities, has not been studied so for in the literature. Unlike gas stations, where the service time is usually within several minutes, the charging period of EV takes several hours, which means that the charging usually occurs at the parking lots. The real-time information of parking lots and gas stations are accessible from the internet in the E-map (e.g., google map) or the online service website of the refueling and parking facilities. Taking the advantage of real-time data from the internet, a practical computation of energy distribution of EV charging load is developed and; therefore, the impact on the power network and the controllability of this additional charging load can be estimated. Load computation of EV charging stations has been reported extensively in the literature. Monte Carlo simulation is used to model the stochastic process of EV charging considering the driving habits and trip distances, etc., and hence to formulate the daily charging power profile of a charging station by summing up the energy consumption of EVs [7,8]. Queuing theory has also been adopted to model the operation of parking lots [9,10] and EV charging demand in [11,12], where the charging processes are represented by an *M*/*M*/∞ queuing system and the arrival rates of EVs are generated as a Poisson process. Vehicle mobility statistics extracted from the 2009 (US) National Highway Travel Survey (NHTS) data is used in [11,12] to define the parameters of the queuing model for the output profiles, for the number of EVs charged and the times EVs arrive and depart from the charging station. However, the stochastic process of EV travel patterns in various countries receives the same input from fixed statistics of a certain area, such as the NHTS dataset. Moreover, the Poisson process in the queuing model assumes EVs have several constant arrival rates and the charging time is typically modeled by Gaussian distribution with given upper and lower limits, which are randomly assigned to each EV. In the proposed queuing theory-based modeling of EV charging stations, the stochastic process of EV arrival and parking can be expressed as a probability function and the parameters can be calculated according to the statistical data accessible from E-map and the service website of facilities.

The modeling and management of EV charging stations are integrated in to the power system analysis and optimal dispatch problem. In the state of art control framework for high-penetration EVs in the power system (i.e., vehicle-to-grid (V2G) technology), the EV aggregator (EVA) is usually introduced as the intermediate control entity between EVs and the distribution system operator (DSO) [13,14]. Since a group of EVs are aggregated at the charging stations, an aggregator is assigned to monitor and control the chargers. The interactive control among EV aggregators and the autonomous control within each aggregator are performed in the typical V2G control framework to relieve the computation load of DSO for the regional distribution network. If properly controlled, the negative impact of the additional charging load can be mitigated and ancillary services can be provided, such as cost minimization, peak shaving, and power quality improvement [12,15]. However, EV parking behavior and energy demand are usually modeled by typical distribution function with the fixed parameters or historical statistics. But the random distribution function and historical statistics of a certain region cannot reflect the diverse characteristics of service facilities for automobiles in other places and; therefore, the results may be very different from the actual condition of the power network with massive EVs.

The control method is the key technique to explore the controllability of charging load for functioning as a virtual energy storage and coordinating with renewable energy sources. There has been a lot of existing research on the optimal dispatching of EV charging power in the power system [16]. From the perspective of an EV aggregator, the location and size of charging aggregations need to be found, the differential evolution and particle swarm optimization (PSO) algorithms are used to solve the optimization problem with the objective of total cost minimization [17]. The location of the charging stations and the optimal number of charging spots were given in the optimal solution [18]. However, the power network has not been taken into account in the model [19]. Other objective functions were investigated for the optimal placement and sizing of an EV aggregator, particularly the profit maximization of an EV aggregator, but some practical issues from the side of users were not described, such as the satisfaction of mobile needs and power quality deterioration caused by recharging massive EVs at peak hours [20]. The optimal dispatch by using EVs as mobile energy storage in the power grid is also studied in the recent literature [21]. In the optimal control problem of EV charging stations modeled in this paper, both EV aggregator configuration at available parking lots and EV charging power regulation should be solved. Therefore, a combined optimization of EV aggregation configuration and operation should be derived based on the modeling of EV charging stations in the optimal dispatch problem.

The main objectives of this paper are to propose a quantitative evaluation of EV charging load and the controllability in the power network, based on the assumption that the fuel consumption at the gas stations is transformed to EV charging facilities covering the same region. The main contributions of the work presented in this paper are as follows:

The operational model of EV charging stations changed from gas stations is presented in Section 2. The optimal configuration and control of charging stations integrated into the optimal dispatch of the power system are explained in detail in Section 3. Section 4 describes the two-stage hybrid optimization algorithm devised for the optimal charging strategy in Section 3. Section 5 presents the analysis and discussion of controllability evaluation for charging facilities in the test 123-bus test network. Section 6 outlines the conclusion drawn based on the numerical results.

#### **2. Modeling of EV Charging Stations Transformed from Gas Stations**

As vehicles served by the gas stations will be replaced by EVs in the future, EV charging facilities will be built to meet this energy demand. The EV charging stations are established with both fast charging and slow charging devices installed at the parking lots. EV charging load changed from gas stations was calculated based on the data extracted from accessible E-map or the online service website of the facilities. For instance, the location and real-time occupancy rate of gas stations and parking lots across the certain area of the town, as shown in Figure 1. The key parameters of the modeling were extracted from the online graphic data to characterize the parking of EVs at different places. The mathematical model of the charging process was formulated considering the energy consumption at gas stations and the EV parking period at parking lots. The weekly occupancy rate of the parking lot, as laid out on the map of Figure 1, shows that each day of the week had a different service pattern. Therefore, a daily occupancy rate was extracted for modeling the operation of EV charging stations, as shown in the right side of the figure. The EV-to-carpark membership function and the queuing theory were also adopted to calculate the parking time and charging service pattern. The temporal and spatial distribution of the EV charging can be calculated according to the online data and the control strategies can be developed based on the adjustable variables, with boundaries defined by modeling the operation of EV charging stations. The workflow of the proposed modeling of the EV charging station and the integrated optimal charging algorithm in the power network is illustrated in Figure 2.

**Figure 1.** Modeling of the electric vehicle (EV) charging facilities changed from gas stations based on online data analysis.

**Figure 2.** Flowchart of the proposed optimal control algorithm for EV charging facilities.

#### *2.1. Oil-To-Electricity Transformation with Membership Function between EVs and Parking Lots*

It was assumed that the fuel consumption at the gas station was transformed to the electricity supplied by EV charging stations during the electrification of transportation system. The energy consumption taken over by the regional EV charging stations can be expressed by:

$$F\_{\varepsilon} = \sum\_{t=1}^{T} \sum\_{j=1}^{G} C\_{f}^{j}(t) = \sum\_{t=1}^{T} \sum\_{k=1}^{M} C\_{\varepsilon}^{k}(t) \tag{1}$$

$$\mathbf{C}\_{\epsilon}^{k}(t) = \sum\_{i=1}^{N^{k}(t)} P\_{i}^{k}(t) \cdot \Delta t \tag{2}$$

where *M* is the number of parking lots, *G* is the number gas of stations covering the same region, *t* denotes *t th* hour, and *T* represents the planned time period, which is a typical day in the simulation, as shown in Figure 1. *j* denotes *j th* gas station, *k* denotes the *kth* parking lot, and *i* is the *i th* EV. *Cf* represents the fuel consumption at the gas station, which is measured by heat based on the quantity (mL) and density (kg/mL) of oil traded at the gas station and fuel heating value (42 MJ/kg), and then converted into the unit of kWh. *Ce* is the equivalent electricity consumption in the unit of kWh. *P*(*t*) is the charging power at each time step, and *Nk*(*t*) is the number of EVs at *kth* parking lot. Because, at each time step, EVs can be at any parking lots in the region, the membership function is defined to represent the possibility and the duration of each EV at the *kth* parking lot. The membership function between EVs to parking lots can be expressed by:

$$\alpha\_{\mathfrak{m}} = \left[ \overleftarrow{a}^1 \overleftarrow{\cdots} \overleftarrow{a}^k \overleftarrow{\cdots} \overleftarrow{a}^{\mathcal{M}} \right] \tag{3}$$

$$\overrightarrow{\boldsymbol{\alpha}}^{k} = \begin{bmatrix} \boldsymbol{\alpha}\_{1}^{k} \cdots \boldsymbol{\alpha}\_{i}^{k} \cdots \boldsymbol{\alpha}\_{N}^{k} \end{bmatrix}^{T} \tag{4}$$

where *Np* is the number of all parked vehicles, α*m* is the matrix of membership degree, and α*k i* is the membership degree of *i th* EV at *kth* parking lot. α*m* subjects to the following conditions:

$$\sum\_{k=1}^{M} \alpha\_i^k = 1 \quad \forall i \in N \tag{5}$$

$$\sum\_{k=1}^{M} N^k(t) = N \quad \forall t \in T \tag{6}$$

$$\sum\_{k=1}^{N} \alpha\_i^k \cdot T = \sum\_{t=1}^{T} N^k(t) \tag{7}$$

These conditions included that the summation of membership degrees for each EV to all the parking lots should be 1. Based on the same assumption that all the vehicles served by the gas stations were changed to parking lots in the same region, the summation of EVs at the *M* parking lots should be equal to the total number of EVs. Furthermore, the total service time of the parking lots can be calculated by the summation of membership degree. It can also be derived from occupancy and capacity of parking lots which can be read from online graphic data as shown in Figure 1. The temporal and spatial distribution of energy consumed at gas stations to the parking lots was calculated by the membership functions in the following section.

#### *2.2. Operation of EV Charging Stations Based on Queuing Theory and Accessible E-Map Database*

The arrival and departure times of EVs at parking lots, which are essential for daily profile of EV load, were derived by the parking data analysis, whereby queuing theory applies. The condition of a typical queuing theory model is that the arrival intensity is constant over the time. But the arrival of EVs for commuting can be concentrated on several hours. Moreover, the service time of parking, which usually takes several hours, is much longer than that of gas station on the time scale of minutes. The steady state of a queuing system cannot be reached in a short time, which is consistent with the real data shown in Figure 1. Thus, in the proposed model, the number of EV was composed of two portions, the steady state and the dynamic state. The steady state was defined as basic EV flow, where the arrival intensity was constant over the time. The dynamic state was defined as categorized EV flow, where EVs were used for a certain purpose and show a fixed pattern, such as household charging for EVs that started charging once coming back from work.

The operation of parking lots was modeled based on queuing theory and the online data. The variation of the parked vehicles can be expressed by:

$$N^k(t) = L^k\_\varepsilon(t) = L^k\_b + N^k\_c(t) \quad \forall k \in M \tag{8}$$

$$
\overline{N\_{\mathcal{L}}(t)} = \overline{L\_{\mathcal{S}}(t)} - \overline{L\_{b}} \tag{9}
$$

where *Lk b* is the expected number of basic EV fleet, which is also the queuing length in the steady state; and *Nk c* is the number of categorized EV fleet, which is varying over the time according to the arrival intensity. These numbers can be read from online service data for each parking lots, as shown in Figure 1. The accessible data is described as *Nc*(*t*), *Ls*(*t*), and *Lb* in Equation (9) for the varying number of EVs parked at each time slot.

According to queuing theory, the arrival of EVs at parking lots is Poisson distribution, denoted by λ*b* and λ*c* for steady and categorized EV fleets, respectively. The queuing length (*Lb*) in the steady state can be calculated by applying the result of queuing theory, where the parking lots can be modeled as a M/M/Z/Z queuing system [7,8].

$$
\lambda\_z = \lambda\_{\mathbb{b}^\star} \text{ z } = 1, 2, \cdots, \; Z^k \tag{10}
$$

$$\mu\_z = z\mu, \; \mathbf{z} = 1, \mathbf{2}, \dots, \; \mathbf{z}^k \tag{11}$$

where *Zk* is the capacity of *kth* parking lot, and *z* is the number of EVs at a parking lot. The service time in this case is the parking period, which complies with negative exponent distribution 1/μ in the queuing model. The service intensity of the parking lot is defined by the queuing theory, and the number of EVs at the parking lot can be written as a function of ρ:

$$
\rho = \frac{\lambda\_b}{Z^k \mu} \tag{12}
$$

$$L\_b^k = \sum\_{z=0}^{z^k} z P\_z = z^k \rho(1 - P\_z) \tag{13}$$

The complete modeling of parking lots is described in Appendix A. The arrival rate of a basic EV fleet can be derived if the queuing length is given:

$$
\lambda\_b = Z^k \mu \rho(\overline{L^k\_b}) \tag{14}
$$

where *Lk b* is the average number of basic EVs at a parking that can be collected from online data.

Similarly, the categorized EV fleet in the queuing model can be quantified by assuming that the arrival rate at each time step is Poisson distribution. The arrival of categorized EVs is concentrated on several hours, which can be easily identified from the occupancy variation of the parking lots. As shown in Figure 1, the rising zone is assumed to be arrival period of categorized EVs as explained in Appendix A. Therefore, the number of arrival EVs *Nck in* can be calculated from *Nk c* (0) and *Nk c* (*t*), as shown in Appendix A, and can also be calculated by:

$$E\left(N\_{\rm ini}^{\rm exk}(t)\right) = \sum\_{y=0}^{\infty} y \frac{\left(\lambda t\right)^k}{k!} e^{-\lambda t} = \lambda\_c^k(t) \cdot \Delta t \tag{15}$$

$$
\lambda\_c^k(t) = \frac{\overline{N\_{\rm in}^k(t)}}{\Delta t} = \frac{\overline{\Delta N\_c^k(t)}}{\Delta t} \tag{16}
$$

where Δ*Nk c* (*t*) is the difference of categorized EV number at *t th* hour, which can be read from the online data shown in Figure 1. The arrival rate at each time slot can be calculated by Equations (14) and (16), and used to determine the parking time period of EVs at parking lots.

The membership degree of each individual EV to parking lot is correlated with the parking period in the queuing theory:

$$
\sigma\_i^k \cdot T = \tau\_i \sim \frac{1}{\mu'} \cdot \frac{1}{\mu} = \overline{\tau\_i} \tag{17}
$$

where τ*i* is the parking period of *t th* EV at one parking lot. The service time complies with negative exponential distribution [7,8]. τ*i* is the average EV parking period that can be collected from online data. As all the other parameters are calculated, the arrival and the time duration of EVs at the parking lot can be completely quantified.

#### *2.3. Parameter Configuration of EV Charging Facilities*

In order to analyze the temporal and spatial distribution of EV charging energy at parking lots covering the same region of gas stations, the number of EVs charged at the parking lot and the energy consumption were calculated according to the fuel consumption at gas stations. Based on the formulation of EV charging facilities installed at parking lots by using membership function, the total number of served EVs can be calculated by:

$$\sum\_{k=1}^{M} \sum\_{i=1}^{N\_{in}^{k}} \alpha\_{i}^{k} \cdot T = \sum\_{i=1}^{N} \sum\_{k=1}^{M} \alpha\_{i}^{k} \cdot T = \sum\_{i=1}^{N} T = NT \tag{18}$$

$$\sum\_{k=1}^{\mathcal{M}} \alpha\_a^k \cdot \mathbf{N}\_{\rm in}^k = \mathbf{N}^P \tag{19}$$

$$a\_i^k = \begin{cases} \begin{array}{ll} \alpha\_a^k & \alpha\_i^k \neq 0\\ 0 & \alpha\_i^k = 0 \end{array} \tag{20}$$

where *NP* is the number of served EVs at all parking lots during *T* hours. The membership degree of EVs at the same parking lot was assumed to be the same in Equation (20), since the customers of a parking lot usually had the same purpose for working in an office building or for staying at home in the residential area. The number of EVs charged at the parking lot can be calculated according to the energy consumption of gas stations:

$$\begin{array}{rcl} F\_{\mathcal{C}} & = \sum\_{k=1}^{M} \sum\_{t=1}^{T} \sum\_{i=1}^{N^k(t)} P\_i^k(t) \cdot \Delta t = \sum\_{k=1}^{M} \sum\_{i=1}^{N^k(t)} \sum\_{t=1}^{T} P\_i^k(t) \cdot \Delta t \\ & = \sum\_{k=1}^{M} \sum\_{i=1}^{N\_{\text{int}}^{k}} EC\_i \Big( SOC\_{out}^k - SOC\_{in}^k \Big) \end{array} \tag{21}$$

where *ECi* is the battery energy capacity of *i th* EV. *SOCk out* is the SOC of EV battery at departure and *SOCk in* is the SOC of EV battery at arrival, respectively. The initial SOC and the target SOC (defined as above 90%) are also random variables that are defined as normal distribution:

$$\text{SOC}\_{\text{in}'}^k \text{SOC}\_{\text{out}}^k \sim \text{N}\left(\mu, \sigma^2\right) \tag{22}$$

The parameters of normal distribution are given for estimating the EV charging load, and the variety of different types of EVs was taken into account:

$$\overline{EC\_i} = \sum\_{r=1}^{R} EC\_r \cdot \eta\_r \tag{23}$$

*ECi* is the average battery capacity. *ECr* is the battery capacity of *r* type EV, and η*r* is the ratio of *r* type EVs to all EVs served at the regional charging stations.

$$F\_{\varepsilon} = \sum\_{k=1}^{M} \sum\_{i=1}^{N\_{\text{in}}^{k}} \overline{EC\_{i}} \cdot \left( \overline{SOC\_{\text{out}}^{k}} - \overline{SOC\_{\text{in}}^{k}} \right) \tag{24}$$

$$\overline{N} = \frac{F\_{\varepsilon}}{\overline{EC\_{\text{i}}} \overline{\{SOC\_{out}^{k}} - \overline{SOC\_{in}^{k}}}} \tag{25}$$

$$
\eta\_{\rm EC} = \frac{\overline{N}}{N^P} \tag{26}
$$

where *N* is the average number of charged EVs, and η*EC* is the ratio of charged EVs to parked EVs.

#### *2.4. Mathematical Formulation of EV Charging Load and Controllable EV Charging Variables*

In order to analyze the impact of EV charging transferred from gas stations and to develop the control strategy for the EV charging process, the calculation of EV charging load and the control variables were derived from the operational model of aforementioned charging facilities incorporated with parking lots.

$$\mathbf{F}\_{\varepsilon}^{\text{sur}} = \sum\_{k=1}^{M} \sum\_{t=1}^{T} \mathbf{C}\_{\varepsilon}^{k}(t) = \sum\_{k=1}^{M} \sum\_{i=1}^{N} \alpha\_{m} \cdot \overleftarrow{P\_{r}} \cdot T \tag{27}$$

$$\overrightarrow{P\_r} = \begin{bmatrix} P\_{r\_1}^1 \cdots \; \; \; \; \; \; \; P\_{r\_2}^k \cdots \; \; \; \; \; P\_r^M \end{bmatrix}^T \tag{28}$$

*Fun e* denotes the energy consumption calculated in an uncontrolled case, and *Pk r* is the rated power of EV charger at *kth* parking lot. The matrix of membership degree of EV to parking lots can also be expressed as α*i* = α1 *i* , α2 *i* , ··· , α*k i* , ··· , α*M i T* .

$$\alpha\_i^k \cdot P\_r^k = \begin{cases} \alpha\_i^k \cdot P\_r^k & \alpha\_i^k \le \frac{T\_i^{kr}}{T} \\ \frac{P\_r^k \cdot T\_r^{kr}}{T} & \alpha\_i^k > \frac{T\_i^{kr}}{T} \end{cases} \tag{29}$$

$$T\_i^{k\tau} = \frac{EC\_{i^\tau}(\text{soc}\_{out}^k - \text{soc}\_{in}^k)}{P\_r^k} \tag{30}$$

*Tkr i* is the charging period of EV to reach target *SOC* at the departure time with rated power.

$$\overrightarrow{P\_{\rm var}^{k}}(t) = \begin{bmatrix} P\_1^k(t), \dots, P\_i^k(t), \dots, P\_N^k(t) \end{bmatrix} \tag{31}$$

*Energies* **2019**, *12*, 1577

$$F\_{\varepsilon}^{\text{cou}} = \sum\_{k=1}^{M} \sum\_{t=1}^{T} \mathbb{C}\_{\varepsilon}^{k}(t) = \sum\_{k=1}^{M} \sum\_{i=1}^{N} \sum\_{t=1}^{T\_{ai}} P\_{i}^{k}(t) \Delta t = \sum\_{k=1}^{M} \sum\_{i=1}^{N} \sum\_{t=1}^{T\_{ai}} P\_{\text{rar}}^{k}(t) \Delta t \tag{32}$$

$$T\_{ai} = a\_i^k \cdot T \tag{33}$$

$$\to \tag{33}$$

*Fcon e* denotes the energy consumption calculated in controlled case, and *Pk* var is the vector of variables (charging power) in a controlled case.

$$\sum\_{t=1}^{T\_{ai}} P\_i^k(t) \cdot \Delta t = E \mathbf{C}\_i \cdot \left( S \mathbf{O} \mathbf{C}\_{out}^k - S \mathbf{O} \mathbf{C}\_{in}^k \right) \tag{34}$$

*T*α*i* is the time period of *i th* EV parked at *kth* parking lot.

#### **3. Coordinated Integration of EV Charging Stations into the Electric Distribution Grid**

As temporal and spatial distribution of EV charging demand and controllable variables are given from the above quantitative analysis of oil-to-electricity transformation for regional automobiles, the charging facilities can be integrated into the electric distribution network. The impact of EV charging load on the power system operation is quantified and the control strategy is developed by adjusting the control variables within the boundaries that are given in the above section. In the vehicle-to-grid operation framework, an EV aggregator is usually defined as the control entity between power grid and EVs, which offers services to aggregate EV charging power (or discharging defined as V2G) from a group of EVs, and acts towards the grid to provide considerable power regulation capacity. In this context, each parking lot is correlated with one EV aggregator, which can perform real-time monitoring and regulation on the charging behavior of each individual EV as long as it is plugged into the power grid. Since the power feedback from EV to grid results in serious battery degradation, which may prevent EV owners from participating in the vehicle-to-grid interaction, only the charging adjustment is taken into account in this paper to provide the optimal charging plan for EVs.

In the impact analysis of EV charging load, EVs were assumed to start charging as soon as they are parked, which is the direct charging used as the uncontrolled case in the simulation. In this case, EVs are charged at rated power of the charging facility or on-board battery charger. The daily charging profile in the impact analysis was defined as the reference to be compared with the controlled case in the simulation. A two-stage optimization was formulated for managing the charging load based on the control variables and the boundaries defined in the above section.

#### *3.1. Objective Function of the Optimal Control Strategy*

The multi-period optimal dispatch with EV charging infrastructure, formulated in considering network constraints, is a mixed integer nonlinear programming (MINLP) problem [20–22]. A two-stage optimization algorithm was derived from the basic formulation of the MINLP optimal charging problem. The optimal configuration of the EV aggregator and the operation of the EV charging infrastructure for regulating the charging rate of individual EVs were modelled in the optimal dispatch problem. The overall objective in this work is to minimize the operating cost. In the objective function, the expenditure on supplying electricity to the load, including charging EVs and the total power loss, were calculated in Equation (35):

$$\min \sum\_{t=0}^{T} \sum\_{i=1}^{n} \rho(t) (P\_{Ld}(t) + P\_{EVAi}(t) + P\_{Ls}(t)) \Delta t \tag{35}$$

where ρ(*t*) is the electricity price at time *t*; *PLd*(*t*) is the net power load of the original distribution system; *PEVAi*(*t*) is the accumulated EV charging power at bus *i*; *PLs*(*t*) is the total power loss at time *t*; and Δ*t* is the time interval, which was defined as one hour in the proposed day-ahead optimal dispatch

of EV charging infrastructure, and it was also consistent with the above charging station model. The load of the EV charging station was calculated by:

$$P\_{EVAi}(t) = \sum\_{m=1}^{Mi} P\_{EVAi}^m(t) \tag{36}$$

*Pm EVAi*(*t*) is the charging rate of individual EVs in the control realm of the EV aggregator installed at bus *i*.

$$P\_{LS}(t) = \sum\_{(i,j)\in B} R\_{ij} \frac{P\_{ij}^2 + Q\_{ij}^2}{|V\_i|^2} \tag{37}$$

*Pij* and *Qij* are the active and reactive power, respectively, transmitted on the branch between bus *i* and *j*; |*Vi*| is the voltage magnitude at bus *i*, and *Rij* is the resistance of branch *ij*.

#### *3.2. Constraints of EV Charging Optimization*

The mathematical model was subject to constraints from each participant in the EV charging optimization, including limits on the reliable operation of power system, EV charging station, and the requirements of EV users. The power flow balance was expressed by:

$$\left(I\_{\vec{i}\vec{j}}(t)\right)^{2} = \left(G\_{\vec{i}\vec{j}}^{2} + B\_{\vec{i}\vec{j}}^{2}\right) \left[\mathcal{U}\_{\vec{i}}^{2}(t) + \mathcal{U}\_{\vec{j}}^{2}(t) - 2\mathcal{U}\_{\vec{i}}(t)\mathcal{U}\_{\vec{j}}(t)\cos\theta\_{\vec{i}\vec{j}}(t)\right] \le I\_{\vec{i}\vec{j}\max}^{2}\tag{38}$$

*Iij*(*t*) is the branch current; *Gij* and *Bij* denote the admittance matrix of the power network; *Ui* and *Uj* are the bus voltage; and *Iij*max is the maximum current of the transmission line.

$$\left(\left(V\_i^{\min}\right)^2 \le V\_i^r(t)^2 + V\_i^{im}(t)^2 \le \left(V\_i^{\max}\right)^2\tag{39}$$

*V*min *i* and *V*max *i* are the lower and upper limits of bus voltage, respectively; *Vr i* (*t*) and *Vim i* (*t*) are the real and imaginary parts of the bus voltage, respectively.

The adjustable EV charging power was limited within the capacity of the charging facility:

$$P\_{EVAi,\min} \le P\_{EVAi}^{\text{m}}(t) \le \min \left\{ P\_{m,\max\prime}, P\_{EVAi,\max} \right\} \tag{40}$$

*PEVAi*,min and *PEVAi*,max are the lower and upper power boundaries of charging facilities, respectively; *Pm*,max is the rated power specified by the EV on-board charger.

The target SOC set by EV users must be achieved till the departure of the parking lots:

$$\text{SOC}\_{\text{tar}}^{\text{un}} \le \sum\_{t=1}^{T} P\_{EVAi}^{\text{m}}(t) \eta\_{clm} / \text{C}\_{bat}^{\text{m}} + \text{SOC}\_{\text{ini}}^{\text{un}} \le 1 \tag{41}$$

where *Cm bat* is the battery capacity of the *mth* EV; *SOCm ini* and *SOCm tar* are the initial and target SOC, respectively; and η*chr* denotes the efficiency of battery charging.

#### **4. Two-Stage Hybrid Optimization Algorithm**

In the proposed EV charging optimization model, the decision variables included the setting of a charging station associated with its EV aggregator and the multi-period charging power control under the given configuration parameters. For the time-series simulation, the increased number of variables made it even harder to solve the MINLP problem. Thus, a two-stage algorithm was formulated to separate the variables into two categories, which were EV aggregator configuration and optimal charging plan for individual EVs. The original problem was converted to a master problem of optimal configuration and a sub-problem of optimal charging strategy, so as to reduce the complexity of the nonlinear constrained optimization modelled in the above section.

#### *4.1. Solution Procedure for Optimal Integration of EV Charging Stations into the Power Grid*

A two-stage hybrid optimization algorithm based on PSO and sequential quadratic programming (SQP) was derived to accelerate the calculation speed, as depicted in Figure 2. As shown in the flowchart, the problem solution of the hybrid optimization algorithm required an iterative process between the master problem of optimal configuration and the sub-problem of optimal charging strategy. In the master problem, the optimal locations of EV charging stations with the power level limits were found by applying PSO algorithm, and the charging plans for each individual EV were embedded as a sub-problem solved in the lower stage. In this way, the aggregated power of EVs connected to the bus of the distribution network was provided so that the PSO algorithm was implemented on mid-voltage level with a limited number of variables. The optimal charging stations were chosen from the candidate locations that were all the parking lots in the region of oil-to-electricity transformation. The capacity of the charging station must be lower than the parking lots:

$$LC\_{EVAi} \in \left\{ LP\_{V1\prime}, \dots, LP\_{Vj\prime}, \dots, LP\_{Vm} \right\} \tag{42}$$

$$CP\_{EVAi} \le (N\_{in}^{ck} + L\_b^k) P\_{EVri} \tag{43}$$

*LCEVAi* denotes the location of charging station at bus *i*, which is selected from the candidate locations of all parking lots, and the bus number is denoted by *LPVj*. *CPEVAi* denotes the capacity of the charging station at bus *i*, and *PEVri* is the rated charging power of charging facility at this charging station.

The capacity of the installed charging facilities was also subject to the operation limits of the power grid considering each EV charging station connected to a bus of the network, as described in the above constraints Equations (38) and (39). All served EVs calculated in Equation (25) were distributed proportionally to each charging station based on its capacity, as expressed by:

$$N\_{EVAi} = \overline{N} \cdot \bigcirc\_{EVAi} / \sum\_{i=1}^{n} \bigcirc\_{EVAi} \tag{44}$$

$$P\_{i, \text{min}} \le P\_{LDi}(t) + \sum\_{m=1}^{N\_{EVAi}} P\_{EVAi}^m(t) \le P\_{i, \text{max}} \tag{45}$$

*NEVAi* is the number of EVs distributed to a charging station according to the proportion of capacities of all charging stations. *PLDi*(*t*) is the base load at bus *i*, *Pi*,min and *Pi*,max are the minimum and maximum power allowed to be connected to bus *i* by the power system, respectively.

The optimal solution of the master problem was provided to the sub-problem as input parameters, including the bus where the EV charging station connected to the network and the charging power limits. The sub-problem without any integer variables can be solved by SQP within the feasible region given in the master problem. The constraints from the power network operation and EV mobile needs were also checked to ensure that the charging plan of each individual EV can fulfill the overall objective.

#### *4.2. Two-Stage Hybrid Optimization Algorithm*

The objective function of the master problem is described in Equation (35), with the location and power capacity of EV charging station in the power network as decision variables and subject to:

$$0 \le P\_{EVAi}(t) \le P\_{EVAi}^{\max}(t) \tag{46}$$

where *P*max *EVAi*(*t*) is the upper limit of the overall charging power, that is the power capacity of each EV aggregator.

$$P\_{EVAi}^{\text{max}}(t) = \sum\_{m=1}^{Mi} P C\_{EVAi}^{\text{m}}(t) \tag{47}$$

*PCm EVAi*(*t*) is the power capacity of individual EV at bus *i*.

$$L\mathcal{L}\_{EVAi} \in \{0, 1\} \tag{48}$$

*LCEVAi* is the location of EV aggregation in the charging station, and in this algorithm is defined as a binary variable.

$$N\_A^{\min} \le \sum\_{i=1}^n L\mathcal{C}\_{EVAi} \le N\_A^{\max} \tag{49}$$

*N*min *A* and *N*max *A* are the minimum and maximum number of EV aggregators planned to be installed in the power system depending on the construction budget of the EV charging infrastructure.

The sub-problem implemented the optimal dispatch of EV charging power under the configuration parameters of the EV aggregator given in the master problem. Accordingly, the sub-problem was formulated as:

$$\min \text{SP}\left(P\_{EVAi}^{\text{pr}}(t)\Big|\_{P\_{EVAi}^{\text{max}}(t), L\subset\_{EVAi}}\right) = \sum\_{t=0}^{T} \sum\_{i=1}^{n} P(t) \left(P\_{Ld}(t) + P\_{EVAi}^{\text{m}}(t) + P\_{Ls}(t)\right) \Delta t \tag{50}$$

where *Pm EVAi*(*t*) *P*max *EVAi*(*t*),*LCEVAi* represents the charging rate of individual EV, which is the decision variable of the sub-problem that is to be optimized under the parameters given in the master problem.

$$\mathbf{f}\_{\rm mf}(L) = \frac{1}{\sqrt{2\pi}\sigma\_{\rm mf}L} \exp\left(\frac{-\left(\ln L - \mu\_{\rm mf}\right)^2}{2\sigma\_{\rm mf}^2}\right) \tag{51}$$

f*m*(*L*) is the probability distribution function of the distance *L* travelled by each EV; σ*m* and μ*m* are the parameters of exponential distribution that is used to simulate the stochastic travelled distance of EV. In this paper, it can be set to ln *L* ∼ *N*(3.46, 0.952).

On the basis of travel distance, the required power recharged of each EV was calculated by:

$$\text{SOC}^{\text{m}}\_{\text{ini}} = \text{SOC}^{\text{m}}\_{\text{tar}} - \frac{\lambda L}{\text{C}^{\text{m}}\_{\text{bat}}} \tag{52}$$

where λ is the energy consumption per unit of distance.

The other commonly used constraints that describe the limits on the charging rate and EV battery SOC boundaries during the charging process under the control of EV aggregator can be found in existing research work [20–23].

#### **5. Simulation Results of the EV Charging Station in Oil-To-Electricity Transformation**

The simulation model built for oil-to-electricity transformation for regional automobiles was mainly composed of EV charging facilities changed from gas stations and the regional distribution network across a small town, as shown in Figures 1 and 3.

#### *5.1. Case Study*

The IEEE 123-bus distribution network with three-phase unbalanced load was used for simulation studies [24], as depicted in Figure 3. For time-series simulation, active and reactive loads were calculated by multiplying the coefficient that reflected the daily load curve, as shown in Figure 4a. The net load profile of the power network was calculated based on the wind power and solar power generation and average daily load of South Australia in [25]. The original load of the test power network was equivalent to the load level indicated by coefficient 0.8. The two-tariff pricing was adopted for daily electric energy price including two main periods. The electricity prices applied to charging stations during the peak hours and the off-peak hours were set to 1.01 and 0.25 CNY/kWh, according to the trial operation of EV charging stations in Shenzhen, Guangdong province of China. It can be seen in Figure 4a that the integration of solar and wind power generation exacerbated the load variation. EV charging load can be controlled for load leveling as shown in Figure 4b, thus ensuring the reliability of the power network. The minimal number and the maximal number of EV charging stations planned to be installed in the test network were set to 5 and 8, considering existing parking lots in the region where automobiles were served by gas stations. The operating cost of the EV aggregator was assumed to be zero in the objective function.

**Figure 3.** EV charging facilities integrated into the 123-bus distribution network.

**Figure 4.** Power load profiles in the distribution network with renewable energy sources: (**a**) Load profiles of 123-bus power network; (**b**) EV charging load with and without control.

Charging stations transferred from two gas stations, as shown in Figure 1, were integrated into the regional distribution network represented by the test 123-node systems shown in Figure 3. Two of the charging stations were assumed to be the parking lots located in the residential area. The characteristics of EV arrival and parking period were identified as the home charging scenario. Two of the charging stations were assumed to be located among office buildings where most of the customers came to work. The last charging station was assumed to be in the commercial area, and composed of EVs parked for shopping and EVs parked overnight. By integrating charging stations into the regional distribution network, the two home charging stations were connected to bus 87 and bus 107 in the 123-node power system. The two office charging stations were connected to bus 31 and bus 39, and the complex charging station was connected to bus 1. The capacity of charging stations, which also means the charging load distribution among the five stations, will be calculated in the proposed EV charging optimization algorithm. In this paper, the acceptable range of bus voltage was set from 0.95 to 1.05 p.u. as the constraints of optimal configuration for EV charging stations. Four different types of EVs were adopted in the case study. As shown in Table 1, the battery capacity and charging rate varied among different EV models. The power electronic interfaced chargers were installed to meet the fast and slow charging needs of EVs, which ranged from 2.2 to 44 kW, and the efficiency of energy conversion process was set to 0.9 in the simulation. As type I and II were designed to be capable of fast charging, these two types of EVs were assigned to office and commercial charging stations. Type III and IV were designed for slow charging and can be connected to standard household outlets. Therefore, these two types of EVs were assigned to home and also the complex charging station where both fast and slow charging facilities were installed. The parameters of EV charging scenario for computing battery energy and charging load are given in Tables 1 and 2.

**Table 1.** Simulation parameters for EV type and EV charging infrastructure.

**Table 2.** Optimal configuration of EV charging stations.

#### *5.2. Results and Discussion*

In the formulation of EV charging stations in oil-to-electricity transformation, charging load distributed at certain charging stations was equivalent to energy consumption for the automobiles served at the gas stations. The optimal configuration of charging stations is given in Table 2. The averages of normal distribution for initial SOC and target SOC at minimum were set to 0.5 and 0.9. The total number of EVs served by the charging stations was calculated to be 420 vehicles, with the assumption that the proportion of different types of EVs at each charging station was the same. For the office charging station at bus 31, the total number of EVs for type I and II was 20. Similarly, 150 EVs of type III and IV were served by the home charging station at bus 107. Four types of EV were served by the complex charging station at bus 1: 40 EVs of type I and II and 100 EVs of type III and IV. The number of EVs at each station were calculated according to the ratio of the capacities given by the proposed optimization algorithm for all EV charging stations. As shown in Table 2 and Figure 3, five optimal locations were selected from the parking lots in the regional power network and the capacities were set to minimize the impact of the additional EV charging load.

The load pattern and voltage deviation of the network with EVs in the uncontrolled condition and charging optimization were compared to evaluate the efficiency of the proposed method, as shown in Figures 4b and 5a. In both cases, the configuration of EV charging stations adopted the results of the optimization algorithm, thus upgrading the performance of the uncontrolled case. Consequently, the uncontrolled charging with the best configuration was picked up from other randomly configured charging stations to analyze the advantage of optimally controlled EV charging load. The charging behavior of EV aggregations was regulated to achieve the valley filling and peak shaving, as seen in Figure 4b. The total power loss was reduced since the power losses at the peak hours were significantly curtailed, as shown in Figure 5b, in which ten lines with largest power flow were selected. The median, minimum, and maximum values of line losses can be seen in the box plot for the ten lines. The bus voltage of the test three-phase unbalanced power network is given in Figure 6, showing the reduction in voltage deviation. For the home charging stations at bus 87 and 107, most EVs were parked overnight and charged for commuting. The comparison of load and voltage profiles indicated that EV charging load is switched to the off-peak period during 24:00 to 6:00 the next day as shown in Figure 7a. Similar results can be found in Figure 7b for the office charging stations at bus 31 and 39, where the two types of EVs were charged to match the network load profile as long as the constraints on the target SOC and the parking time can be satisfied. Noted that the two-tariff pricing directed the EV charging load to off-peak hours with a lower electricity price to cut down the cost, and the reduced power loss further minimized the objective function.

**Figure 5.** Power losses of distribution network with EV charging load: (**a**) Daily profile of total losses with and without EV charging control; (**b**) line losses of 123-bus power network at the peak hour.

**Figure 6.** Bus voltage of the 123-bus test system with and without EV charging control.

**Figure 7.** Load and voltage profiles at bus 31 phase-C and bus 107 phase-B: (**a**) Load profiles at different buses; (**b**) voltage profiles at different buses with and without EV charging control.

To calculate the charging plan for each EV, the capacity of the charging station fixed by the optimal configuration restricted the maximal charging power at each bus. The overall charging power during the parking period was allocated to each individual EV by the EV aggregator at a charging station. The charging plans for EV type I in uncontrolled and controlled cases are compared in Figure 8a, and the corresponding SOC curves for all EVs of type I are depicted in Figure 8b. The charging plans of EVs connected at different buses are shown in Figure 9, in which the charging profiles of two types of EVs served at the charging station connected to node-31 phase-C and node-107 phase-B are illustrated. The slow charging of type III and IV needed several hours to reach the target SOC, while the fast charging for type I and II can be finished within an hour if the parking period of the vehicle was limited. The charging rate of each individual EV at both fast charging and slow charging facilities was coordinated to enhance the economic and reliable operation of power system, without breaking the boundaries on service capacity of facilities and EV mobile needs described in Sections 2 and 3. It can be seen that EVs of different types reached the desired SOC before departure even though the charging power profile varied according to the solution given by optimal charging algorithm.

**Figure 8.** Charging rate and state of charge (SOC) profiles of EV type I at fast charging facilities with and without control: (**a**) EV charging curves (Tesla Model X); (**b**) SOC curves during the parking period.

**Figure 9.** Charging profiles of EVs connected to facilities at different buses: (**a**) Charging profiles of EV type-I and type-II parked at bus 31; (**b**) charging profiles of EV type-III and type-IV parked at bus 107.

#### **6. Conclusions**

To evaluate the impact of charging stations changed from gas stations in the oil-to-electricity transformation for automobiles, a mapping relation model between gas stations and charging stations serving vehicles in the same region was established by using membership degree function of EV and queuing model-based charging stations. A novel feature of the proposed EV charging stations was that the real-time data of gas stations and parking lots serving the regional automobiles was used to determine the stochastic EV charging time and charging demand distribution in the power network. Therefore, the controllability of EV charging load can be quantified by integrating the EV charging process with parameters and control variables in the optimal dispatch of the distribution network. With the defined boundaries and power system constraints, a combined optimization of configuration and operation of EV charging stations was developed for reducing the operating cost and risks of power system. A two-stage hybrid optimization algorithm composed of master problem and sub-problem was proposed to reduce the complexity of the combined optimization problem. The optimal setting of charging stations and the optimal charging created for each EV can improve the economic and reliable operation of the power system with flexible EV charging load. The proposed method can serve as an effective tool for evaluating the increased EV charging load and the controllability of charging stations in the oil-to-electricity transformation for automobiles across the town.

**Author Contributions:** Methodology and writing—original draft preparation, S.G.; writing—review and editing, J.W.; project administration, B.X.

**Funding:** This research was partially funded by the National Key R&D Program of China (2016YFB0900400), and the National Nature Science Foundation of China (51707130).

**Acknowledgments:** This research was partially supported by the International Clean Energy Talent Programme (iCET). The three-month training in Sweden improved the research work in the paper and contributed to the data collection for modeling the regional gas station and parking lots from their online service websites and electronic map.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A Detailed Modeling of Parking Lots Based on Queuing Theory**

This section presents the modeling of parking lots as a M/M/Z/Z queuing system [7,8]. The variation of the parked vehicles given the data accessible online can be expressed by:

$$N^k(t) = L^k\_\mathfrak{s}(t) = L^k\_\mathfrak{b} + N^k\_\mathfrak{c}(t) \qquad \forall k \in \mathcal{M} \tag{A1}$$

$$
\overline{N\_{\mathbb{C}}(t)} = \overline{L\_{\mathbb{S}}(t)} - \overline{L\_{b}} \tag{A2}
$$

$$N\_b^k(t+1) = N\_b^k(t) + \Delta N\_b^k(t+1), \qquad \Delta N\_b(t+1) = 0 \qquad \forall t \in T \tag{A3}$$

$$N\_{\varepsilon}^{k}(t+1) = N\_{\varepsilon}^{k}(t) + \Delta N\_{\varepsilon}^{k}(t+1), \qquad N\_{\varepsilon}(0) = 0 \qquad \forall t \in T \tag{A4}$$

$$
\Delta \mathcal{N}\_b^k(t) = \mathcal{N}\_{\rm in}^{\rm bk}(t) - \mathcal{N}\_{\rm out}^{\rm bk}(t) \tag{A5}
$$

$$
\Delta \mathcal{N}\_c^k(t) = \mathcal{N}\_{in}^{ck}(t) - \mathcal{N}\_{out}^{ck}(t) \tag{A6}
$$

where *Lk b* is the expected number of the basic EV fleet, which is also the queuing length in the steady state, and *Nk c* is the number of categorized EV fleet, which is varying over the time according to the arrival intensity.

The total number of both EV fleet can be read from the online data, but the arrival number of EVs at each time needs to be calculated. According to queuing theory, the arrival of EVs at parking lots is Poisson distribution:

$$N\_{\rm in}^{\rm c}(t) \sim \lambda\_{\rm c} \tag{A7}$$

$$N\_{\rm in}^b(t) \sim \lambda\_b \tag{A8}$$

where λ*b* and λ*c* are the average arrival rate of the basic and categorized EVs respectively. The queuing length *Lb* in the steady state can be calculated by applying the result of queuing theory, the parameters of a M/M/Z/Z queuing system are adopted for parking lots:

$$
\lambda\_z = \lambda\_b, \text{ z } = 1, 2, \dots, \text{ } \text{ } Z^k \tag{A9}
$$

$$\mu\_z = z\mu,\;\mathbf{z} = 1,2,\dots,\;Z^k\tag{A10}$$

where *Zk* is the capacity of *kth* parking lot, and *z* is the number of EVs at a parking lot. The service time in this case is the parking period, which complies with negative exponent distribution 1 μ in the queuing model.

$$P\_0 = \left(\sum\_{i=0}^{z^k} \frac{\left(z^k e\right)^i}{i!} \right)^{-1} \qquad \text{z} = 1, 2, \cdots, \; Z^k \tag{A11}$$

$$P\_z = \frac{\left(z^k \rho\right)^z}{z!} P\_0 \tag{A12}$$

$$
\rho = \frac{\lambda\_b}{Z^k \mu} \tag{A13}
$$

*P*0 denotes the probability of no EV at parking, and *Pz* denotes the probability of *z* EVs at parking. ρ is the service intensity of the parking lot defined by the queuing theory. Consequently, the number of EVs at the parking lot can be written as a function of ρ, and the arrival rate of basic EV fleet can be derived if the queuing length is given.

$$L\_b^k = \sum\_{z=0}^{z^k} z P\_z = z^k \rho (1 - P\_{z^k}) \tag{A14}$$

$$L\_b^k = \overline{L\_b^k} \tag{A15}$$

where *Lk b* is the average number of basic EVs at parking that can be collected from online data.

Similarly, the categorized EV fleet in the queuing model can be quantified by assuming that the arrival rate at each time step is Poisson distribution. The arrival of categorized EVs is concentrated on several hours which can b easily identified from the occupancy variation of the parking lots. As shown in Figure 1, the rising zone is assumed to be arrival period of categorized EV. Therefore, the arrival EV number *Nck in* can be calculated from *Nk c* (0) and *Nk c* (*t*) as shown in Equations (A2)–(A6), and can also be calculated by:

$$N\_{in}^{ck}(t) = N^k(t) - L\_b^k \tag{A16}$$

$$N\_{\rm in}^{\rm ck} = \sum\_{t=1}^{T\_{\rm c}} N\_{\rm in}^{\rm ck}(t), \qquad N\_{\rm in}^{\rm ck}(t) \sim \lambda\_{\rm c}^{k}(t) \tag{A17}$$

$$\overline{N\_{in}^{ck}(t)} = \lambda\_c^k(t) \cdot \Delta t \qquad t \in [0, T\_c] \tag{A18}$$

where *Nck in*(*t*) is the average number of categorized EVs arriving at *t th* hour, which can be collected from online data as shown in Figure 1. λ*k c* (*t*) is the Poisson distribution of categorized EVs arriving at *t th* hour. *Tc* is the time period of categorized EV arrival. Following the queuing model of each parking lot, the arrival rate of categorized EV fleet at each time step can be derived by:

$$P(N\_{\rm in}^{\rm ck}(t) = \mathbf{y}) = \frac{(\lambda t)}{y!} e^{-\lambda t} \qquad \lambda = \lambda\_{\mathbf{c}}^k \qquad \mathbf{y} = 0, 1, 2, \cdots \tag{A19}$$

$$N\_{in}^{k}(t) = N\_{in}^{ck}(t) + N\_{in}^{pk}(t) \tag{A20}$$

*Nk in* is total number of EVs parked at *kth* parking lot within *Tc* hours, and *Nk in*(*t*) is the number of EVs arriving at *kth* parking lot during each time interval. As described above, it is composed of two portions. *Nbk in* (*t*) denotes the number of basic EVs arriving at *kth* parking lot at *t th* hour, and *Nck in*(*t*) denotes the number of categorized EVs arriving at *kth* parking lot at *t th* hour. Based on the queuing theory, the number of EVs arriving at each time step is given by:

$$E\left(N\_{\rm in}^{\rm rk}(t)\right) = \sum\_{y=0}^{\infty} y \frac{\left(\lambda t\right)^{k}}{k!} e^{-\lambda t} = \lambda\_{c}^{k}(t) \cdot \Delta t \tag{A21}$$

$$
\lambda\_c^k(t) = \frac{\overline{N\_{\rm in}^k(t)}}{\Delta t} = \frac{\overline{\Delta N\_c^k(t)}}{\Delta t} \tag{A22}
$$

where Δ*Nk c* (*t*) is the difference of categorized EV number at *t th* hour that can be read from the online data shown in Figure 1.

The membership degree of each individual EV to parking lot is correlated with the parking period in the queuing theory:

$$
\alpha\_i^k \cdot T = \tau\_{i\prime} \qquad \tau\_i \sim \frac{1}{\mu} \tag{A23}
$$

where τ*i* is the parking period of *i th* EV at one parking lot. The service time complies with negative exponential distribution:

$$P(\tau\_i \le t) = \begin{cases} 1 - e^{-\mu t} & t \ge 0 \\ 0 & t < 0 \end{cases} \tag{A24}$$

$$E(\tau\_i) = \int\_0^{+\infty} \frac{+\infty}{0} \,\mu + e^{-\mu t} dt = \frac{1}{\mu} \tag{A25}$$

$$\frac{1}{\mu} = \overline{a\_i^k} \cdot T = \overline{\tau\_i} \tag{A26}$$

where τ*i* is the average EV parking period that can be collected from online data. As all the other parameters are calculated, the arrival and the duration of EVs at the parking lot can be completely quantified.

#### **References**

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


**Author details**

Patricia Acosta-Vargas1 , Janio Jadán-Guerrero2 , Cesar Guevara2 , Sandra Sanchez-Gordon3 \* and Tania Calle-Jimenez3

1 Intelligent and Interactive Systems Lab, Universidad de Las Américas, Quito, Ecuador

2 Universidad Tecnológica Indoamérica, Quito, Ecuador

3 Department of Informatics and Computer Science, Escuela Politécnica Nacional, Quito, Ecuador

\*Address all correspondence to<EMAIL_ADDRESS>
© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**93**

*Technical Contributions to the Quality of Telerehabilitation Platforms: Case Study—ePHoRt…*

[9] International Organization for Standardization. ISO 9241- 11—Ergonomics of human system

Vol. 208. 2015. pp. 221-225

Boundaries; 2013

[10] Kushniruk A, Monkman H, Tuden D, Bellwood P, Borycki E. Integrating heuristic evaluation with cognitive walkthrough: Development of a hybrid usability inspection method. In: Studies in Health Technology and Informatics.

[11] Xiao L, Yan X, Emery A. Design and evaluation of web interfaces for informal care providers in senior monitoring. In: Proceedings of the 76th ASIS&T Annual Meeting: Beyond the Cloud: Rethinking Information

[12] Gulati A, Sanjay K. Critical analysis on usability evaluation techniques. International Journal of Engineering Science and Technology (IJEST). 2012;**4**

[13] Lim C, Hae-Deok S, Lee Y. Improving the usability of the user interface for a digital textbook platform

for elementary-school students. Educational Technology Research Development. 2012;**60**:159-173

[14] Quinones D, Rusu C. How to develop usability heuristics: A systematic literature review. Computer Standards and Interfaces. 2017;**53**:89-122

[15] Pilco H, Sanchez-Gordon S, Calle-Jimenez T, Rybarczyk Y, Jadán J, Villarreal S, et al. Analysis and improvement of the usability of a telerehabilitation platform for hip surgery patients. In: Advances in Intelligent Systems and Computing. Vol. 781. Cham: Springer; 2019. pp. 197-209

[16] Acosta-Vargas P, Acosta T, Lujan-Mora S. Challenges to assess accessibility in higher education websites: A comparative study of

interaction. 2018

*DOI: http://dx.doi.org/10.5772/intechopen.83686*

[1] Rautaray S, Agrawal A. Vision based hand gesture recognition for human computer interaction: A survey. Artificial Intelligence Review.

[2] Ichim A, Bouaziz S, Pauly M. Dynamic 3D avatar creation from hand-held video input. ACM Transactions on Graphics.

[3] Pavone E, Tieri G, Rizza G, Tidoni E, Grisoni L, Aglioti S. Embodying others in immersive virtual reality: Electrocortical signatures of monitoring the errors in the actions of an avatar seen from a first-person perspective. The Journal of Neuroscience.

[4] Belal S, Taktak A, Nevill A, Spencer S. A fuzzy system for detecting distorted plethysmogram pulses in neonates and paediatric patients. Physiological Measurement. 2001;**22**(2):397-412

[5] Guevara C, Jadán-Guerrero J, Rybarczyk Y, Acosta-Vargas P, Esparza W, González M, et al. A real-time algorithm for movement assessment using fuzzy logic of hip arthroplasty patients. In: Advances in Intelligent Systems and Computing. Vol. 781. Cham: Springer; 2019. pp. 265-273

[6] Han J, Shao L, Xu D, Shotton J. Enhanced computer vision with Microsoft Kinect sensor: A review. IEEE Transactions on Cybernetics.

[7] Sastry S. Intelligent Control. CRC

[8] Suganthi L, Iniyan S, Samuel A. Applications of fuzzy logic in renewable energy systems—A review. Renewable and Sustainable Energy

Reviews. 2015;**48**:585-607

2013;**43**(5):1318-1334

Press; 2018

2015;**34**(4):45:1-45:14

2016;**36**(2):268-279

2015;**43**(1):1-54

**References**

*Technical Contributions to the Quality of Telerehabilitation Platforms: Case Study—ePHoRt… DOI: http://dx.doi.org/10.5772/intechopen.83686*


8. Startability

The Kantrowitz criterion for intake starting [26] says that the normal shock, in front of an intake duct, will move downstream and out of the duct if the flow at the exit of the duct is not choked—the duct flow will start. This criterion applies to the normal shock at the entry of the duct as well as at any other position in the duct. On a wavecatcher intake, Figures 18 and 19, overboard flow spillage will allow the shock to move downstream, over the external/open portion of the intake, until it reaches the V-notch at the beginning of the internal flow. This is made possible by

The Busemann Air Intake for Hypersonic Speeds DOI: http://dx.doi.org/10.5772/intechopen.82736

an effectively large flow area, on the moving, post-shock side, allowing overboard flow spillage. If the Kantrowitz condition, for the shock at the point of inflection, at the V-notch, is satisfied the shock will continue moving downstream, out through the internal flow section, and the intake will start.

The entry area of the internal flow, Af, is defined as the conical surface at the angular position where the surface is inflected because the flow is normal to the conical surface at the inflection and a stationary, conical, normal shock is compatible with the flow there. The size of this area is available from a Busemann intake calculation. This area is needed for application of the Kantrowitz starting criterion.

## 8.1 Startability of the weak shock Busemann

Note that, for given M2 and δ23, Eq. (9), gives two solutions for θ23, for a weak and a strong shock. This leads to the possibility of generating two different Busemann intakes, the weak shock version would have supersonic and the strong shock would have a subsonic exit flow. Because of its supersonic exit flow the weak shock intake is better suited for scramjet application. However, at contractions to be useful for scramjets, the Busemann intake with a weak shock does not start spontaneously or if it does start then it does so for intakes with an insufficient amount of contraction.

The determination of startability for a wavecatcher Busemann intake is as follows. At first we examine the startable weak-shock Busemann flow to show that it does not provide sufficient compression:

Many such calculations, starting from weak shock waves, (os) in Figure 3, were performed with the outcomes plotted on a graph of area ratio, A3/A1 vs. entry Mach number, M1, in Figure 24. Each result is shown as a dot that is colored green if the totally internal flow Busemann intake duct starts, green or yellow if the wavecatcher Busemann intake module starts (as determined in d) above) and red if there is no start. Curves of the "startability index," S = (A1 – Ai)/(AK – Ai), measure the location of a dot on the overall area ratio scale where S = 0 on the isentrope and S = 1 on the Kantrowitz criterion. Intermediate, fractional values, are on curves between these limits. The curve for S = 0.6 seems to well represent the startability limit for wavecatcher Busemann intake designs based on the weak shock condition. As seen from the figure the wavecatcher design lowers the startable area ratio from about 0.6 to 0.4. This is still not good enough. For good engine performance, it is

Intake performance, both efficiency and capability, are not affected much by considerable amounts of either truncation or stunting. This is due to the fact that the high-loss leading edge boundary layer flow is not eliminated but merely moved downstream. Also there is some increase in inviscid flow losses from the finite angle leading edge from truncation. However, an estimated 15–30% weight saving

is available through wall materials elimination resulting from intake fore-

Variation of intake compression and contraction as caused by truncation and stunting.

Hypersonic Vehicles - Past, Present and Future Developments

were found in [10].

Figure 23.

8. Startability

104

7.4 Leading edge blunting

are formed, keeping the exit flow uniform.

shortening. It appears that the significant advantage of truncation and stunting is not to intake performance but to the saving of structural weight. Similar results

Busemann flow has no deflection at the leading edge so that the leading edge tends to be sharp and thin. Such leading edges are difficult to cool at hypersonic speeds. Transpiration cooling is made possible by a slight rounding of the leading edge. Rounding or blunting affects both the viscous as well as the inviscid flow in the intake [10]. The strong bow shock causes a hot entropy layer to overlay the boundary layer and cause it to thicken. The same shock focuses on the symmetry axis producing a Mach reflection at the center line. It was shown in [10] that a 1 mm diameter leading edge on a 500 mm diameter Busemann intake, flying at Mach 10 and 30 km altitude, is optimal in reducing the viscous and inviscid losses. It seems that the combination of blunting and stunting should be such that the conical shock is kept incident on the Busemann surface corner, so that no reflected shock waves

The Kantrowitz criterion for intake starting [26] says that the normal shock, in front of an intake duct, will move downstream and out of the duct if the flow at the exit of the duct is not choked—the duct flow will start. This criterion applies to the normal shock at the entry of the duct as well as at any other position in the duct. On a wavecatcher intake, Figures 18 and 19, overboard flow spillage will allow the shock to move downstream, over the external/open portion of the intake, until it reaches the V-notch at the beginning of the internal flow. This is made possible by

back-pressure is lowered to be replaced by a weak shock structure. This now opens a possibility. We calculate a Busemann intake module, with a strong shock, such that it starts spontaneously, at a high overall contraction ratio, and then reduce the back-pressure to remove the strong shock and obtain a supersonic exit flow, with a weak shock. This yields an exit flow which is suitable for scramjet operation. In doing this, we note that the flow, from the freestream to the location of the strong shock has not changed as we switch to the weak shock, so that the intake remains on the strong-shock design flow up to the corner while being started. Also, the amount of internal contraction remains the same and we could really start the intake with the weak shock structure in the first place. The strong Busemann shape is really a design tool which leads to a modified Busemann flow but with a started intake of high compression and efficiency having a supersonic exit Mach number-an intake with a high overall contraction but with a low internal, startable contraction. The supersonic exit flow is no longer conical although its axial symmetry is preserved. In comparison to the weak shock case (Figure 24) there is a considerable enhancement of startability in the enlarged yellow domain so that starting can be expected near the S = 0.1 curve, which is acceptable for scramjet applications.

The reason for this improvement in startability stems from the fact that, for the strong shock option, the angular distance between the strong shock (at the corner) and the inflection cone is small so that Af and A3 areas are close in size, i.e., there is not much internal contraction. This makes it easy for the conical normal shock, at the inflection location, to be swallowed. The strong shock design's aim is to produce a wavecatcher Busemann intake with a high overall contraction ratio and a low internal contraction that starts spontaneously. The calculated shape is compatible with a normal conical shock positioned at the inflection location (angle) and we select the intakes, with internal contractions, that permit the shock to be swallowed. Flow downstream of the inflection shock is subsonic. The strong shock is present only fleetingly during flow starting. After swallowing a weak shock system appears in the exit with supersonic flow downstream. This flow is suitable for a scramjet combustor. Within the calculated strong shock contour the supersonic weak shock flow is no longer conically symmetric and has to be examined via CFD and experiment. The end result is a wavecatcher Busemann intake with a high overall contraction (S 0.1) but with a low internal, self-starting contraction and hence a

Figure 26 is a schlieren picture of four quarter-circle Busemann intake modules each with flows started by the unsteady gun tunnel starting flow [4]. Arrow points

Impulsively started 4-module Busemann intake in the gun tunnel at Mach 8.33 (model is similar to

startable intake and supersonic exit flow.

The Busemann Air Intake for Hypersonic Speeds DOI: http://dx.doi.org/10.5772/intechopen.82736

Figure 26.

Figure 19).

107

#### Figure 24.

Busemann intake startability based on weak shock design: green—full Busemann start; green and yellow wavecatcher start with spillage; and red—no start.

desirable to operate intakes near the curve S = 0.1 in Figure 24. Startable weak shock Busemann designs do not produce this amount of contraction—the prospect of obtaining startable full or wavecatcher intakes, with sufficient compression, from weak shock Busemann flows is disappointing.

#### 8.2 Startability of the strong shock Busemann

Startability calculations for the full and wavecatcher Busemann intakes were done also by starting the integration of Eqs. (5) and (6) from strong shock waves, following the (a–e) steps above.

The outcomes are plotted in Figure 25. Each result is shown as a dot that is coloured green if the totally internal flow Busemann intake duct starts, green or yellow if the wavecatcher Busemann intake module starts (as determined in d) above) and red if there is no start.

The strong shock version has high overall contraction but low internal contraction so that it will self-start at overall contractions useful to the scramjet as a wavecatcher, but the strong shock and subsonic exit flow are not useful to the scramjet engine. In a wavecatcher module the about-to-start strong shock will be held in place by the appropriate back-pressure. It will move downstream if the

#### Figure 25.

Busemann intake startability regimes based on strong shock design. green—full Busemann start; green and yellow—module start with spillage; and red—no start.