Volume 16 / Issue 2 SEPTEMBER 2023 www.fe.um.si/en/jet.html 2 JET JET 3 VOLUME 16 / Issue 2 Revija Journal of Energy Technology (JET) je indeksirana v bazah INSPEC© in Proquest’s Technology Research Database. The Journal of Energy Technology (JET) is indexed and abstracted in database INSPEC© and Proquest’s Technology Research Database. 4 JET JOURNAL OF ENERGY TECHNOLOGY Ustanovitelj / FOUNDER Fakulteta za energetiko, UNIVERZA V MARIBORU / FACULTY OF ENERGY TECHNOLOGY, UNIVERSITY OF MARIBOR Izdajatelj / PUBLISHER Fakulteta za energetiko, UNIVERZA V MARIBORU / FACULTY OF ENERGY TECHNOLOGY, UNIVERSITY OF MARIBOR Glavni in odgovorni urednik / EDITOR-IN-CHIEF Jurij AVSEC Souredniki / CO-EDITORS Bruno CVIKL Miralem HADŽISELIMOVIĆ Gorazd HREN Zdravko PRAUNSEIS Sebastijan SEME Bojan ŠTUMBERGER Janez USENIK Peter VIRTIČ Ivan ŽAGAR Uredniško izdajateljski svet / PUBLISHING & EDITORIAL COUNCIL Dr. Anton BERGANT, Litostroj Power d.d., Slovenia Prof. dr. Marinko BARUKČIĆ, Josip Juraj Strossmayer University of Osijek, Croatia Prof. dr. Goga CVETKOVSKI, Ss. Cyril and Methodius University in Skopje, Macedonia Prof. dr. Nenad CVETKOVIĆ, University of Nis, Serbia Prof. ddr. Denis ĐONLAGIĆ, University of Maribor, Slovenia JET 5 Doc. dr. Brigita FERČEC, University of Maribor, Slovenia Prof. dr. Željko HEDERIĆ, Josip Juraj Strossmayer University of Osijek, Croatia Prof. dr. Marko JESENIK, University of Maribor, Slovenia Prof. dr. Ivan Aleksander KODELI, Jožef Stefan Institute, Slovenia Prof. dr. Rebeka KOVAČIČ LUKMAN, University of Maribor, Slovenia Prof. dr. Milan MARČIČ, University of Maribor, Slovenia Prof. dr. Igor MEDVED, Slovak University of Technology in Bratislava, Slovakia Prof. dr. Matej MENCINGER, University of Maribor, Slovenia Prof. dr. Greg NATERER, Memorial University of Newfoundland, Canada Prof. dr. Enrico NOBILE, University of Trieste, Italia Prof. dr. Urška LAVRENČIČ ŠTANGAR, University of Ljubljana, Slovenia Izr. prof. dr. Luka SNOJ, Jožef Stefan Institute, Slovenia Prof. Simon ŠPACAPAN, University of Maribor, Slovenia Prof. dr. Gorazd ŠTUMBERGER, University of Maribor, Slovenia Prof. dr. Anton TRNIK, Constantine the Philosopher University in Nitra, Slovakia Prof. dr. Zdravko VIRAG, University of Zagreb, Croatia Prof. dr. Mykhailo ZAGIRNYAK, Kremenchuk Mykhailo Ostrohradskyi National University, Ukraine Prof. dr. Marija ŽIVIĆ, Josip Juraj Strossmayer University of Osijek, Croatia 6 JET Tehnični urednik / TECHNICAL EDITOR Sonja Novak Tehnična podpora / TECHNICAL SUPPORT Tamara BREČKO BOGOVČIČ Izhajanje revije / PUBLISHING Revija izhaja štirikrat letno v nakladi 100 izvodov. Članki so dostopni na spletni strani revije www.fe.um.si/si/jet.html / The journal is published four times a year. Articles are available at the journal’s home page - www.fe.um.si/en/jet.html. Cena posameznega izvoda revije (brez DDV) / Price per issue (VAT not included in price): 50,00 EUR. Informacije o naročninah / Subscription information: http://www.fe.um.si/en/jet/subscriptions.html Lektoriranje / LANGUAGE EDITING TAIA INT d.o.o. Oblikovanje in tisk / DESIGN AND PRINT Tiskarna Saje d.o.o. Naslovna fotografija / COVER PHOTOGRAPH Jurij AVSEC Oblikovanje znaka revije / JOURNAL AND LOGO DESIGN Andrej PREDIN Ustanovni urednik / FOUNDING EDITOR Andrej PREDIN Izdajanje revije JET finančno podpira Javna agencija za raziskovalno dejavnost Republike Slovenije iz sredstev državnega proračuna iz naslova razpisa za sofinanciranje domačih znanstvenih periodičnih publikacij / The Journal of Energy Technology is co-financed by the Slovenian Research Agency. JET 7 Spoštovani bralci revije Journal of energy technology (JET) Zgodovina energetike najverjetneje sega tisočletja nazaj, v čas, ko je človek spoznal moč ognja. Ta je nato s pomočjo ognja, biomase in nekaterih drugih energentov ustvaril kvalitetna orodja, energijo pa je uporabljal za pripravo hrane in ogrevanje prostorov. Tisoče let pred našim štetjem so določene civilizacije že posedovale znanje o namakanju polj, v Mezopotamiji in starem Egiptu pa so se celo zavedali pomena vodnih ekosistemov. Tako so civilizacije ob tokovih Evfrata, Tigrisa in Nila pričele z razvojem poljedelstva in živinoreje. V času industrijske revolucije ob koncu 18. stoletja so se zgodile velike spremembe na področju porabe energije, saj se je takrat začela intenzivna izraba premoga kot energenta, intenzivna izraba nafte pa se je začela nekoliko kasneje, v začetku 19. stoletja. To je bilo obdobje mnogih odkritij na področju energetike, razvoja elektrarn oz. pretvorbe energije v električno energijo, tehnologije motorjev, materialov, … Intenzivna izraba energentov je žal povzročila velike ekološke probleme, hkrati pa se je še povečalo število svetovnega prebivalstva. Na začetku drugega tisočletja se je pričela nova revolucija v energetiki, ki ji pravimo internetna revolucija. Ob razvoju internetnih tehnologij smo priča sodobnim energetskim procesom z visoko stopnjo učinkovitosti. V tem desetletju pa smo že v drugi revoluciji, ki jo zaznamuje razvoj umetne inteligence. Uporaba umetne inteligence v energetiki je sicer še v povojih, toda zaradi zmožnosti obdelave velike količine podatkov se trenutno že dogaja nepredviden razvoj. Upajmo le, da bodo vse te novosti koristile tako človeku kot naravi. Jurij AVSEC odgovorni urednik revije JET 8 JET Dear Readers of the Journal of Energy Technology (JET) The history of energy can be traced back millennia, when humankind first came to realise the power of fire. With this, humans created quality tools to make their lives easier, with the help of fire, their surrounding nature, and some other energy sources, using the energy to cook food and heat rooms. Even in the millennia before Common Era, there was already substantial knowledge in certain societies regarding land irrigation, and the importance of water ecosystems was already a staple of Mesopotamia and ancient Egypt. Thus, communities along the Euphrates, Tigris, and Nile rivers embarked upon the cultivation of crops and the domestication of animals. A seismic shift in energy practices unfolded during the Industrial Revolution at the close of the 18th century, marked by the widespread adoption of coal as a predominant energy source. Subsequently, the early 19th century witnessed the emergence of intensive oil usage, accompanied by a multitude of breakthroughs in the field of energy science. This era witnessed the establishment of power plants, and innovations in electricity, engines, technology, and materials. Unfortunately, with the population of the planet increasing more than ever before, the intensive use of energy products has caused major ecological problems. At the beginning of the second millennium, a new revolution in energy began, otherwise known as the Internet Revolution. With the development of internet technologies, we have witnessed major modern energy processes with high levels of efficiency. Meanwhile, in the present decade, we are witnessing a new revolution, marked by the development of artificial intelligence. The use of artificial intelligence in the energy industry is still in its infancy, but with its ability to process a large amount of data, the unsuspected development that we are now witnessing is finally occurring. Let us only hope that all these innovations will benefit both humankind and nature. Jurij AVSEC Editor-in-chief of JET JET 9 Table of Contents Kazalo Using household flexibility as a part of demand side management programs Uporaba modela fleksibilnosti gospodinjstva kot dela programov upravljanja na strani povpraševanja Katerina Bilbiloska, Aleksandra Krkoleva Mateska, Petar Krstevski . . . . . . . . . . . . . . . . . . . . . 11 Optimal dispersed generation placement in radial distribution systems using a pattern of load flow procedure – clustering optimisation Optimalna postavitev razpršenih virov energije v radialnih distribucijskih omrežjih po vzorcu določanja pretoka moči – optimizacija grozdenja Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski . . . . . . . . . . . . . . . . . . . . 23 Visualization of 4d thermal maps 4d vizualizacija temperaturnega polja Gorazd Hren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Seawater as a reagent in the flue gas desulphurisation process Morska voda kot reagent v procesu reažvepljevanja dimnih plinov Martin Bricl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Instructions for authors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 10 JET JET Volume 16 (2023) p.p. 11-22 Issue 2, 2023 Type of article: 1.01 http://www.fe.um.si/si/jet.htm USING HOUSEHOLD FLEXIBILITY AS A PART OF DEMAND-SIDE MANAGEMENT PROGRAMS UPORABA MODELA FLEKSIBILNOSTI GOSPODINJSTVA KOT DELA PROGRAMOV UPRAVLJANJA NA STRANI POVPRAŠEVANJA Katerina Bilbiloska1, Aleksandra Krkoleva Mateska R,1Petar Krstevski1 Keywords: demand side management, linear optimisation Abstract This paper presents how customer flexibility can be used in the framework of a Demand Side Management (DSM) program. It focuses on household consumers, and aims to show how slight behavioural changes can lead to optimised use of household appliances and consequently, reduced energy costs. The paper presents a solution that would enable household customers to practice their flexibility. The solution is based on linear optimisation technique (LP) which is used to redistribute the use of domestic appliances throughout the day. The results show that implementation of the proposed solution can lead to reduction of hourly electricity use and reduce the overall electricity costs of the household. The optimised use of appliances reduces the household peak load more than 40% compared to the baseline scenario for all analysed cases. The optimisation implemented on specific households allows for reducing energy costs on average by 30%. R11R Corresponding author: Prof., Aleksandra Krkoleva Mateska, Ss. Cyril and Methodius University in Skopje, Faculty of Electrical Engineering and Information Technologies, Rugjer Boskovic No.18, 1000 Skopje, North Macedonia, Tel.: +389 2 3099 122, E-mail address: krkoleva@feit.ukim.edu.mk 1 Ss. Cyril and Methodius University in Skopje, Faculty of Electrical Engineering and Information Technologies, Rugjer Boskovic No.18, 1000 Skopje, North Macedonia JET 11 Katerina Bilbiloska, Aleksandra Krkoleva Mateska, Petar Krstevski JET Volume 16 (2023) p.p. Issue 2, 2023 Povzetek V tem članeku predstavimo, kako lahko uporabimo prilagodljivost strank v okviru programa za upravljanje povpraševanja (DSM). Osredotočimo se na gospodinjske odjemalce, saj želimo pokazati, kako lahko majhne spremembe v vedenju vodijo do optimizirane uporabe gospodinjskih aparatov in posledično do znižanja stroškov energije. V prispevku predstavimo rešitev, ki bi gospodinjskim odjemalcem omogočila, da prakticirajo svojo fleksibilnost. Rešitev temelji na tehniki linearne optimizacije (LP), ki se uporablja za prerazporeditev uporabe gospodinjskih aparatov čez dan. Rezultati kažejo, da lahko izvajanje predlagane rešitve prispeva k zmanjšanju urne porabe električne energije in tako zmanjša skupne stroške električne energije v gospodinjstvu. Optimizirana uporaba aparatov zmanjša konično obremenitev gospodinjstva za več kot 40 % v primerjavi z osnovnim scenarijem za vse analizirane primere. Optimizacija, izvedena na posameznih gospodinjstvih, omogoča znižanje stroškov energije v povprečju za 30 %. 1 INTRODUCTION The effects of global warming and climate change are already observable. Statistics related to sea level show its constant increase [1] caused by the melting of glaciers and ice sheets as well as the thermal expansion of warm seawater. The Earth is becoming warmer due to greenhouse gases emissions in the atmosphere and the uncontrolled use of fossil fuels. To neutralise the negative environmental impacts, governments around the world are adopting strategies [2] and measures [3] to develop and use low-carbon technologies, and to increase the share of renewable energy sources (RES) in an attempt to secure the path towards decarbonisation. RES, as one of the most common solutions for decarbonisation, is related to additional challenges - the operation of power systems with high-RES penetration requires dealing with various technical challenges caused by intermittent nature of these resources. Power generation technologies based on fossil fuels are used to provide for the base load, but some of them (mainly natural gas) are commonly used to provide the required flexibility for the system. In the near future, with decreased use of fossil fuel-based technologies, RES technologies should also compensate for system flexibility. Traditionally, electricity generation follows consumption. This approach can be revised considering that the load can have certain flexibility and thus, the consumption can change to fit the available generation. This actually defines the idea behind the DSM concept. [4] The demand side flexibility is based on the use of load and storage through the DSM concept, combined with innovative solutions. [5] In general, DSM is defined as measures, activities and programs that establish cooperation between electricity suppliers and consumers by encouraging a change in energy use which differs from their usual habits to achieve economic benefit to both parties. Both residential, [6] especially consumers with smart homes [7], [8] and industrial consumers, as discussed in [9], [10], [11] may benefit from the implementation of this concept. DSM strives to equalise the load curve, i.e. to reduce the occurrence of peaks and the duration of peak periods. That could be achieved through several approaches: demand response (DR), energy efficiency and creating opportunities for energy storage. [4] This paper first analyses DSM programs and then presents a specific application of DSM for households, based on a simple and effective approach that optimises household consumption and redistributes the use of various electrical appliances given certain constraints. The aim is to reduce the peak of the load curve and thus, contribute to a more efficient operation of the 12 JET Using Household Flexibility as a Part of Demand Side Management Programs grid and reduction of the electricity bill for the customer. The financial aspect, reduction of the electricity bill is also described through one simple example. 1.1 DSM Programs The DSM concept includes measures, activities and programs that establish cooperation between suppliers and consumers of electricity. It encourages a change in habits related to energy usage in order to achieve economic benefit for both parties. From the point of view of energy suppliers and operators of transmission and distribution systems, DSM strives to equalise the load curve, i.e. to reduce the occurrence of peaks and duration of periods with peaks. The DSM concept assumes that consumers’ comfort is not compromised [12] and that their normal economic operation is maintained. [13], [14] As mentioned above, depending on the specific way of implementation of DSM, different subcategories could be defined: DR, energy efficiency and energy storage. DR is one of the most used concepts in the DSM framework. It defines change in energy use by consumers compared to their normal consumption, as a response to external signal. DR is an effective approach to use consumers’ flexibility to change the inelastic consumption into elastic, and at the same time, provide benefits for the consumer through predetermined incentives. [15] In this context, the consumers are usually classified in separate groups as industry, commercial consumers and households. [15] For example, industrial consumers can participate in DR programs if their production process is not disturbed, while for commercial consumers the heating and cooling systems or the lighting can be managed in a manner that allows response to price or other signals. There is great potential for the application of DR in the residential sector, but this largely depends on the implementation of smart devices that can automatically respond to changes without the need for manual adjustment. [16] Depending on the way of implementation, DSM programs can be divided into two categories, explicit (incentive) programs and implicit (price-based) programs. Explicit programs aim to change electricity usage by consumers in response to a request or signal from a supplier or other entities, in return for receiving payment for service or other benefits. [17] These programs allow end users to participate in the wholesale, capacity and balancing markets in accordance with the network rules. They can participate in these markets directly or be represented through a specific entity aggregator, created for that purpose. Namely, the aggregators represent a number of customers and use their flexibility to allow them to participate in explicit DR programs. [18] The balancing services become increasingly interesting for explicit demand side flexibility, as discussed in [19]. Furthermore, introducing electric vehicles as providers of demand side flexibility increases the potential benefits of explicit DSM. [20] A recent study [21] showed that industrial and aggregated DSM providers can also participate in the congestion management markets, which are under development across Europe. As a general categorisation, explicit DSM programs encompass direct load control, curtailable service, demand bidding, emergency demand response, capacity market and ancillary market programs. The implicit programs are based on dynamic definition of the price, i.e. the price of electricity changes during various periods of the day. The consumption should respond to the changes in the electricity prices, which should result in lower consumption during the periods of high prices. Compared to explicit programs, the implicit programs do not foresee direct consumer participation in the electricity markets. However, the price-induced customer behaviour impacts electricity markets, as discussed in [22]. The programs are typically offered by the suppliers and allow the customers to change their consumption in relation to the offered tariffs. In this category, time of use tariff, critical load pricing, extreme day price JET 13 Katerina Bilbiloska, Aleksandra Krkoleva Mateska, Petar Krstevski JET Volume 16 (2023) p.p. Issue 2, 2023 programs, extreme day price with peak (critical) load and real time price (variable price) are included. The design of tariffs requires considering a number of factors to provide adequate results in cost reduction, as discussed in [23] and [24]. The practical implementation of these programs depends on several factors including appropriate regulatory framework, developed electricity markets that foster the participation of both load and generation for providing bids and offers in all market timeframes, especially related to reserve capacity and balancing energy. [25] In addition to this, application of appropriate information and communication technologies (ICT) is considered an enabler for the practical implementation of DSM programs. The ICT solutions should be applied both at supplier/aggregator level and customer level. It is important to use cost-effective solutions that would be affordable and simple to use, especially when these programs should be applied at small customers and households. [26] 2 DSM APPLICATION FOR HOUSEHOLDS BASED ON LINEAR OPTIMISATION The integration of DSM in the households is a complex problem primarily due to the different devices used, as well as the specific habits of the users. The main goal when implementing DSM in a household is not to force change the habits of the users, but to motivate their engagement, considering their comfort constraints. Application of DSM for households could be integrated through simple and effective approach that optimises household consumption and redistributes the use of various electrical appliances given certain constraints. The aim is to reduce the peak of the load curve and thus, contribute to a more efficient operation of the grid and reduction of the electricity bill for the customer. The solution presented in this paper distributes the use of certain devices within a 24-hour period, aiming to minimise the household load at each hour of the observed period. 2.1 Method The method used in this paper is based on LP. The LP optimisation technique is implemented to distribute the use of household electric appliances within specific time periods. The LP process considers the user-defined constraints related to various categories of devices. The basic idea of the developed solution is to determine in which period of the day a specific appliance should be used, with aim to reduce electricity costs, to reshape the load curve and to provide certain service for distribution system operator. If LP is used to find an optimal use schedule of household appliances, then that would mean to determine household appliance use during the day in order to minimise the load in each hour of the day. The LP optimisation is based on the Eq. (1.1) - (1.5), [27]: 14 JET Using Household Flexibility as a Part of Demand Side Management Programs For the purpose of this research, the household appliances are classified in three categories depending on their usual period of use, i.e. 1) appliances with constant time of use; 2) appliances with variable time of use; and 3) devices with variable energy and time of use, as described in detail in [28]. The first category covers appliances for which a change of use is not preferred and includes electric stoves, TVs, refrigerators with freezers etc. The second category encompasses appliances that could be used in various periods of the day like water heater, air conditioner, electric vehicle battery charger, etc. The third category includes appliances with shiftable time of use and energy use that changes over time. It includes devices such as washing machines, dishwashers and clothes drying machines. Eq. (1.3) - (1.5) introduce the user-defined constraints and are related to the three categories of devices mentioned above. The optimisation process consists of two parts: the first part, which optimises the use of the appliances from the second category (with variable time of use) based on the described LP optimisation and the second part, which aims to find the best period for the use of the devices from the third category, while achieving the least peak durations. To solve the optimisation problem, a computational program was developed using standard MATLAB package. The input data includes total number of devices, number of devices per category, installed power of each device, and time of use each of the devices. The input file contains a matrix with 24 rows and eight columns. The rows refer to each hour of the day, while the columns represent the devices being considered in the household. Therefore, the number of columns may change for different households according to the number of observed appliances. The elements of the matrix show the state of the device, i.e., whether the device is turned on or not during the considered period. To simplify the optimisation and to reach the optimal solution faster, the elements of this matrix are 1 and 0 – where 1 indicates the devise is ‘on’ and 0 the device is ‘off’. In addition to this, information about the categories of the devices is obtained with the aim to show how many devices belong to each of the categories that were previously described. In a separate matrix column, the data for total engagement of a separate device during the day, installed power of each device as well as the minimum and maximum power of the device are provided, the latter referring devices with variable time and power of use are entered. After reading the input data and placing it in appropriate variables, an optimisation function is called. The optimisation consists of two parts. The first time the program optimises the operation of devices with variable time of use and it is assumed that they can be turned on at any time of the day. This separation is made due to the specificity of devices with variable energy and time of use because their operation should not be interrupted. The optimisation approach of these devices is equal to shortest path algorithms whereby all possible combinations are considered and the combination that causes the minimum peak load is selected as a desirable solution. JET 15 Katerina Bilbiloska, Aleksandra Krkoleva Mateska, Petar Krstevski 2.2 JET Volume 16 (2023) p.p. Issue 2, 2023 Optimisation results In this paper, the operation of the devices in a four-member household is optimised through the application of a DR program. The input data used for the devices are indicative, based on [28]. It is assumed that the air conditioning unit can be switched on in an economic mode (predefined temperature) during the day to keep a steady temperature in the household for inhabitants returning before 17 hours (usually school children older than 10 years) and it is placed in the second category of appliances. It could easily be switched to another category if the household members have other habits. This is valid for the other appliances as well – their use will depend on the habits of the household members and thus vary from one household to another. The results are presented in Fig. 1 to Fig. 4. Fig. 1 shows the baseline scenario, without any optimisation. Fig. 2 presents the case after the optimisation is applied, showing that in the optimised solution the peak in the household consumption is significantly decreased. To broaden the analyses, additional cases have been considered. The sensitivity analyses applied on the proposed solution show how changes in user defined constraints influence the optimised solution. Fig. 3 shows the load curve if the devices from the second category are switched on only after 5 p.m., while Fig. 4 shows the load curve if a constraint on the number of devices switched on per hour is introduced (in this case, four devices may be simultaneously switched on). Further sensitivity analyses are presented in [29]. It is worth mentioning that this method places the appliances in the solution once a pre-defined condition is satisfied, which explains the occurrence of the peak in early morning hours or after 5 p.m. in the examples below. Figure 1: Baseline scenario - without optimisation 16 JET Using Household Flexibility as a Part of Demand Side Management Programs Figure 2: Optimised use of the considered devices Figure 3: Category 2 devices are switched on after 5 p.m. Figure 4: Load curve with a limited number of switched devices per hour JET 17 Katerina Bilbiloska, Aleksandra Krkoleva Mateska, Petar Krstevski 2.3 JET Volume 16 (2023) p.p. Issue 2, 2023 Financial Impact of Implementation of DSM programs The practical implementation of a DSM program is plausible if both sides – suppliers (aggregators) and customers are expecting to obtain benefits. It is reasonable that the financial compensation (or financial benefit) is the main motivation for both parties to participate in DSM programs. To valorise the expected benefits for the customers, the examples described in the previous subsection are compared to the baseline scenario, considering existing prices of electricity for each case. The performed calculations show the amount of money spent by the household to cover the electricity costs per day if the household appliances are used as in the baseline scenario and as in the scenarios depicted in Fig. 2 to Fig. 4 (the optimised distribution of use of the devices). The calculations are performed considering the current electricity prices for households in the Republic of North Macedonia [30], i.e. low tariff of 1.3183 denars/kWh for the period between 10 pm and 7 am and first block tariff of 4,7257 denars/kWh for the rest of the day (for monthly usage of up to 200 kWh). The prices are for energy only and costs for transmission and distribution of electricity are not included. The results of the calculations are presented in Table 1 below. The period T1 represents the period of the day between 7 am and 10 pm and the period T2 represents the period between 10 pm and 7am. Although the tariff system in place considers 4 tariff blocks, the calculations are done considering the tariff for the first block only. It is a reasonable approach as the optimisation should maximise the use of shiftable appliances in the off-peak period, leading to a situation where the customer monthly consumption would not exceed the limit of the first block. This is valid in the cases where electricity is not the major source for household heating. Furthermore, the calculations are only done once per day, to relate to the simulation results of the previous section. Table 1: Costs for energy for the baseline and optimised scenarios Time of day Energy (kWh) MKD Total (MKD) Total (EUR) Baseline scenario (Fig. 1) Optimised scenario (Fig. 2) Optimised scenario (Fig. 3) Optimised scenario (Fig. 4) T1 58 274.1 T1 31.6 149.3 T1 33.2 156.9 T1 35.4 167.3 T2 1.8 2.37 276.5 4.5 T2 28.2 37.18 186.5 3 193.3 3.14 T2 27.6 36.4 T2 25.4 33.5 200.8 3.3 As can be noted from the results presented in Table 1, the redistribution of the use of the appliances to periods with lower electricity price could significantly decrease the electricity bill. The optimisation presented in Fig. 2 results with the least energy costs, which is achieved in accordance with the customer’s constrains, i.e. without considerable changes of their habits. Relieving customer’s constraints would lead to achieving higher savings. Fig. 5 gives a graphical presentation of the possible reduction of costs relative to the baseline scenario. 18 JET Using Household Flexibility as a Part of Demand Side Management Programs 100% 90% 80% 70% 60% 50% 40% 30% 20% 10% 0% Baseline scenario (Fig. 1) Optimized scenario Optimized scenario Optimized scenario (Fig. 2) (Fig. 3) (Fig. 4) Figure 5: Reduction of costs for the analysed scenarios The described example presents the possible household savings on electricity if a specific DR program is in place, with pre-defined time periods and tariffs per periods. In such case, the implementation of the described LP method would allow the customer to use the shiftable household appliances in periods with low tariff following a schedule derived from the LP solution. Further analyses are required to valorise the benefits for customers in the case of dynamic prices. 3 CONCLUSIONS This paper shows how the implementation of simple optimisation solutions can ease the implementation of DSM programs in households. By introducing user-defined constraints, solutions can be tailored towards the consumer, but without the need to implement complex algorithms and advanced technologies. Both the system operators and the consumers benefit from the optimisation of the household flexibility. The presented analyses confirm a decrease in the household peaks, something which benefits the distribution system operation. Furthermore, the distribution of the appliances following the described LP optimisation will lead to decrease of household electricity costs. The savings will depend on the constraints introduced by the customers and their willingness to slightly change their living habits. It is worth mentioning that the implementation of such solutions requires certain behavioural changes, something which are often hard to change. Reduction of electricity costs may be perceived as a strong motivation factor, especially in periods of crisis and high prices, especially in developing economies. However, to maintain implementation, it is also important to raise public awareness on the costs related to electricity production, transmission and distribution and to promote energy efficiency as a priority measure for both public and private sectors. It could be assumed that the benefits of such programs would be higher if households became prosumers and even had their own storage systems. Further analyses are required to confirm and quantify this assumption. JET 19 Katerina Bilbiloska, Aleksandra Krkoleva Mateska, Petar Krstevski JET Volume 16 (2023) p.p. Issue 2, 2023 References [1] Climate gov: https://www.climate.gov/news-features/understanding-climate/climatechange-global-sea-level (29.05.2023) [2] European Commission: Communication from the Commission, The European Green Deal, COM(2019) 640 final, Brussels, 11.12.2019 [3] European Commission: Clean Energy for all Europeans package https://energy.ec.europa. eu/topics/energy-strategy/clean-energy-all-europeans-package_en, (10.04.2022) [4] P. Warren: A review of demand-side management policy in the UK, Renewable and Sustainable Energy Reviews, Vol 29.(C), pp. 941-951, 2014 [5] C. D. Korkas, M. Terzopoulos, C. Tsaknakis, E. B. Kosmatopoulos: Nearly optimal demand side management for energy, thermal, EV and storage loads: An Approximate Dynamic Programming approach for smarter buildings, Energy and Buildings, Vol. 225, 111676, 2022 V. Stavrakas, A. Flamos, A modular high-resolution demand-side management model to quantify benefits of demand-flexibility in the residential sector, Energy Conversion and Management, Vol. 205, 112339, 2020 [6] E. Sarker, P. Halder, M. Seyedmahmoudian, E. Jamei, B. Horan, S. Mekhilef, A. Stojcevski, Progress on the demand side management in smart grid and optimization approaches, International Journal of Energy Research, Vol. 45, pp. 36– 64, 2021 [7] M. Afzal, Q. Huang, W. Amin, K. Umer, A. Raza and M. Naeem, Blockchain Enabled Distributed Demand Side Management in Community Energy System With Smart Homes, IEEE Access, vol. 8, pp. 37428-37439, 2020 [8] A. Dzyuba, I. Solovyeva, Price-based Demand-side Management Model for Industrial and Large Electricity Consumers, International Journal of Energy Economics and Policy, 2020, 10(4), 135-149 [9] H. Dunkelberg, T. Heidrich, T. Weiß and J. Hesselbach, Energy demand flexibilization of industrial consumers, Journal of Simulation, Vol. 14:1., 53-63, 2020 [10] C. Dang, J. Zhang, C. -P. Kwong and L. Li, Demand Side Load Management for Big Industrial Energy Users Under Blockchain-Based Peer-to-Peer Electricity Market, IEEE Transactions on Smart Grid, vol. 10, no. 6, pp. 6426-6435, 2019 [11] S. Nojavan, K. Zare: Demand Response Application in Smart Grids – Operation Issues, Vol. 2, Springer Cham, pp. 75-117, 2020 [12] M. H. Shoreh, P. Siano, M. Shafie – khah, V. Loia: A survey of industrial applications of demand response, Electric Power Systems Research, Vol.141, pp. 31-40, 2016 [13] J. Ponoćko, J.V. Milanović, A. Krkoleva Mateska, P. Krstevski, S. Borozan: Existing Approaches to Wide-scale DSM Deployment to Facilitate Transmission Network Flexibility - Results of the Survey in South-East Europe, Proceedings, 2019 IEEE PES Innovative Smart Grid Technologies Europe, 2019 20 JET Using Household Flexibility as a Part of Demand Side Management Programs [14] J. Aghaei, M.-I. Alizadeh: Demand response in smart electricity grids equipped with renewable energy sources: A review, Renewable and Sustainable Energy Reviews, Vol. 18, pp. 64-72, 2013 [15] C. Goldman, M. Reid, R. Levy, A. Silverstein: National Action Plan for Energy Efficiency -Coordination of Energy Efficiency and Demand Response, United States Environmental Protection Agency, 2010 [16] A. Abiri-Jahromi, N. Dhaliwal, F. Bouffard: Demand Response in Smart Grids, Integration of Demand Response into the Electricity Chain - Challenges, Opportunities and Smart Grid Solutions, (ISTE Ltd, John Wiley & Sons), pp. 1-9, 2015 [17] P. Sianon: Demand response and smart grids—A survey, Renewable and Sustainable Energy Reviews, vol. 30, pp.461-478, 2014 [18] T. Freire-Barcelo, F. Martin-Martinez, Alvaro Snachez-Miralles: A literature review of Explicit Demand Flexibility providing energy services, Electric Power Systems Research, vol. 209, 107953, 2022 [19] J. Einolander, R. Lahdelma: Explicit demand response potential in electric vehicle charging networks: Event-based simulation based on the multivariate copula procedure, Energy, vol. 256, 124656, 2022 [20] J. Ponoćko, A. Krkoleva Mateska, P. Krstevski: Cross-border DSM as a complement to storage and RES in congestion management markets, International Journal of Electrical Power & Energy Systems, Vol. 148, 108917, 2023 [21] B. Vatandoust, B. B. Zad, F. Vallée, J. -F. Toubeau, K. Bruninx: Integrated Forecasting and Scheduling of Implicit Demand Response in Balancing Markets Using Inverse Optimization, 2023 19th International Conference on the European Energy Market, 2023 [22] J. Freier, V. von Loessi: Dynamic electricity tariffs: Designing reasonable pricing schemes for private households, Energy Economics, vol. 112, 106146, 2022 [23] B. Guo, M. Weeks: Dynamic tariffs, demand response, and regulation in retail electricity markets, Energy Economics, vol. 106, 105774, 2022 [24] Smart Energy Demand Coalition: Explicit Demand Response in Europe - Mapping the Market 2017, (SEDC), April 2017 [25] S. Panda, S. Mohanty, P. Kumar Rout, B. Kumar Sahu, M. Bajaj, H. M. Zawbaa, S. Kamel: Residential Demand Side Management model, optimization and future perspective: A review, Energy Reports, Vol. 8, pp. 3727-3766, 2022 [26] Z. Zhu, J. Tang, S. Lambotharan, W. Chin, Z. Fan: An integer linear programming based optimization for home demand-side management in smart grid, Proceedings of 2012 IEEE PES Innovative Smart Grid Technologies Europe, 2012 [27] K. Bilbiloska: Utilization of consumers’ flexibility in power systems by implementation of demand-side management, Master Thesis, Ss. Cyril and Methodius University, Faculty of EE and IT, Skopje, pp. 38-40, 2021 JET 21 Katerina Bilbiloska, Aleksandra Krkoleva Mateska, Petar Krstevski JET Volume 16 (2023) p.p. Issue 2, 2023 [28] K. Bilbiloska, A. Krkoleva Mateska, P. Krstevski: Optimization of Customer Flexibility within Implicit Demand Side Management Programs, Proceedings of 18th International Conference on the European Energy Market (EEM), Ljubljana, September, 2022 [29] Energy and Water Services Regulatory Commission of the Republic of North Macedonia: Electricity prices https://www.erc.org.mk/page_en.aspx?id=287 (20.04.2023) 22 JET JET Volume 16 (2023) p.p. 23-46 Issue 2, 2023 Type of article: 1.01 http://www.fe.um.si/si/jet.htm OPTIMAL DISPERSED GENERATION PLACEMENT IN RADIAL DISTRIBUTION SYSTEMS USING A PATTERN OF LOAD FLOW PROCEDURE – CLUSTERING OPTIMISATION OPTIMALNA POSTAVITEV RAZPRŠENIH VIROV ENERGIJE V RADIALNIH DISTRIBUCIJSKIH OMREŽJIH PO VZORCU DOLOČANJA PRETOKA MOČI – OPTIMIZACIJA GROZDENJA Angela Popova1, Jovica VuletićR,2Jordančo Angelov2,3Mirko Todorovski 24 Keywords: Distributed generation, Power injection placement, Clustering optimisation, Distribution system, Losses R Corresponding author: Associate Professor, Jovica Vuletić, University of Ss. Cyril and Methodius – Faculty of Electrical Engineering and Information Technologies, Power Systems Department, Rugjer Boskovic 18, P.O. Box 574 1000 Skopje, R. N. Macedonia, Tel.: +3892 3099 125, E-mail address: jovicav@pees-feit.edu.mk 1 Technical University of Munich, Moosacher Str. 81, 80809 Munich, Germany, ge75jit@mytum.de, 2 University of Ss. Cyril and Methodius – Faculty of Electrical Engineering and Information Technologies, Power Systems Department, Rugjer Boskovic 18, P.O. Box 574 1000 Skopje, R. N. Macedonia, Tel.: +3892 3099 125, E-mail address: jordanco@pees-feit.edu.mk, mirko@pees-feit.edu.mk JET 23 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 Abstract Nowadays, distributed generation plays a vital role in distribution systems. It makes an indisputable contribution towards power loss minimisation and voltage profile improvement. To maximise the benefits of distributed generators, their location and size is of crucial importance. This paper describes the use of clustering optimisation (CO) as a highly effective method for determining the optimal placement and sizing of distributed generators. Assessment of the effectiveness from the clustering optimisation method is achieved by its demonstration on a 69-bus and 119-bus distribution systems. Furthermore, the results obtained from implementation of the proposed approach are compared to results from recent studies including analytical approaches, heuristic, and meta-heuristic, as well as mathematical programming algorithms. It is concluded that clustering optimisation is a simple and efficient method in terms of optimal allocation and sizing of distributed generators. It outweighs other approaches in terms of simplicity, results, and computation time. Povzetek Porazdeljena proizvodnja energije dandanes igra ključno vlogo v distribucijskih sistemih, saj nesporno prispeva k zmanjšanju izgub moči in izboljšanju napetostnega profila. Da bi povečali prednosti porazdeljenih virov, sta njihova lokacija in velikost ključnega pomena. Ta članek opisuje uporabo optimizacije združevanja v grozde (CO) kot zelo učinkovito metodo za določanje optimalne postavitve in dimenzioniranja porazdeljenih virov. Oceno učinkovitosti metode optimizacije grozdenja dosežemo z njeno demonstracijo na distribucijskih sistemih z 69 in 119 zbiralkami. Poleg tega rezultate, pridobljene z izvajanjem predlaganega pristopa, primerjamo z rezultati iz nedavnih študij, vključno z analitičnimi pristopi, hevristiko in meta-hevristiko ter algoritmi matematičnega programiranja. Na ta način ugotovimo, da je optimizacija združevanja v grozde preprosta in učinkovita metoda v smislu optimalne alokacije in dimenzioniranja porazdeljenih virov ter odtehta druge pristope v preprostosti, rezultatih in času izračuna. 1 INTRODUCTION Increasing climate change, dwindling resources and greenhouse gas emissions have led to an increase in the exploitation of renewable energy sources for electricity production. As renewables are scattered around the country, their potential can mostly be tapped through integration to the distribution system as a form of distributed generation (DG). In recent years, the share of DGs in power systems has been significantly increasing. Distributed generation can be defined as any electricity generating technology, installed by the utility system or at the site of a utility customer, connected at the distribution system level of the electric grid. [1] Therefore, DG integration undoubtedly affects the power system, especially at the distribution level. Properly located and sized units have the potential to reduce total power losses in the system, improve the voltage profile, enhance reliability, enable greater available capacity for power transmission and reduce equipment stress. [2] On the other hand, if not optimally placed and sized, the installation of DGs in the grid could cause an increase of system losses, crippling of the voltage conditions, voltage flicker and an increased level of harmonics, all leading to considerably greater costs. [3] Hence, the use of an efficient optimisation method for sizing and placement is of immense importance to fully access the benefits of DGs. [4-7] 24 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation Numerous technologies can be considered as DGs. Photovoltaics and solar-thermal units, wind turbines, small hydro plants, geothermal units, all types of fuel cells and battery storage technologies can be categorised as DGs. [8] However, from a power system analysis point of view, any type of DG technology, depending on its specifics, can be presented as an active, reactive, or apparent power injection (PIN). Over the past few years, a variety of different approaches and methods have been developed in relation to the issue of optimal sizing and allocation of PINs. They can generally be classified as analytical methods [8-14], heuristic and meta-heuristic [15-28], hybrid methods [29-37] and mathematical programming algorithms [38-41]. All approaches have unique advantages, but also have various limitations and setbacks to some extent. Analytical techniques perform mathematical analysis of power distribution systems, resulting in a set of numerical equations that are used to formulate an objective function [13,14]. They have some remarkable features such as simplicity, easy implementation, and short computation time [7,10]. These techniques are suitable for variety of objectives, but they also have certain limitations. For instance, the authors in [8] present an analytical approach for determining the optimal location and size of a DG. The proposed algorithm consists of formulating a loss sensitivity factor based on equivalent current injection without requiring the Jacobian matrix, the admittance (or its inverse) matrix and manages to give competent results. However, the optimisation procedure only determines the location and size of a single DG unit that only injects active power. Another analytical method is presented in [11] which also considers only active power injection for the purpose of determining maximum power loss reduction. Heuristic methods are mainly based on engineering experience and knowledge acquired through research as done in [16]. These methods aim to explore the search space in a particularly convenient manner. Heuristic methods are generally problem dependent. Their greatest setback, however, is that they always give a possible solution which is not necessarily optimal. For example, the heuristic method based on Uniform Voltage Distribution Algorithm presented in [18] despite being robust and fast, fails to give the optimal allocation of reactive power injections. Meta-heuristic approaches (also known as population-based optimisation methods) are widely used because of their problem-independence and ability to provide competent results. However, most meta-heuristic algorithms are only approximation algorithms as they cannot always find the global optimal solution. Furthermore, they require tuning of a great number of parameters as well as long computation time due to their iteration-based nature. Teaching learning-based optimisation - TLBO [42] algorithm overcomes the problem of defining algorithm specific parameters. The optimality of the results and convergence properties of the TLBO algorithm have been improved in [42] by introducing a new quasi-oppositional TLBO - QOTLBO method. The QOTLBO algorithm; however, is parameter-dependent and has the tendency to get trapped in a local optimum. These shortcomings have been overcome with the comprehensive TLBO - CTLBO method presented in [43]. Nevertheless, all these methods can be classified as meta-heuristic approaches, hence retaining the advantages and disadvantages from this group of methods. Combining methods from different groups, i.e., hybrid methods can contribute to overcoming some of the shortcomings from the individual methods. Therefore, the authors in [29] present a hybrid analytical and meta-heuristic method for optimal placement and sizing of multiple DGs of JET 25 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 different types. DG sizes are determined by analytical method, whereas their locations by a PSO algorithm. The results obtained from this method are superior compared to results attained from PSO and analytical methods. However, due to the iteration-based nature of the PSO method, the algorithm can encounter difficulties in terms of convergence speed and accuracy by increasing the number of variables and implementation on large scale networks. LSF-BFOA method is presented in [31] and used for simultaneous placement of DG and DSTATCOM. Location of the DG and DSTATCOM is determined by using a loss sensitivity factor - LSF and the optimal size is then determined by using a Bacterial Foraging Optimisation Algorithm - BFOA. Despite offering superior results, the LSF-BFOA method is characterised with a relatively long computation time. Something similar is done in [35] using loss sensitivity factors - LSF and then gradually placing DGs in a set determined from the sensitivity analysis. Table 1: Overview and number of methods the proposed approach is compared with Citation Approach AM [8] analytical № of methods citation outperforms 2 IA [10] analytical 2 AA [13] analytical EA [14] NH [15] BPSO-SLFA Method Test case Year IEEE-69 2009 IEEE-69 2013 2 IEEE-69 2015 analytical 5 IEEE-69 2016 heuristic 5 IEEE-69&119 2019 [26] heuristic 3 IEEE-69 2020 TLBO [42] meta-heuristic 4 IEEE-69 2014 QOTLBO [42] meta-heuristic 4 IEEE-69&119 2014 IWO [27] meta-heuristic 5 IEEE-69 2016 CTLBO [43] meta-heuristic 2 IEEE-69&119 2018 DE-TLCHS [25] meta-heuristic 1 IEEE-69 2019 MRFO [28] meta-heuristic 8 IEEE-69 2020 HPSO [29] hybrid 2 IEEE-69 2014 LSF-BFOA [31] hybrid 3 IEEE-119 2016 MFO [37] hybrid / IEEE-69 2017 BIBC [36] hybrid / IEEE-69 2018 QODELFA [32] hybrid 6 IEEE-69&119 2019 LSF-SSA [33] hybrid 2 IEEE-69 2019 DGPI [35] hybrid 6 IEEE-69 2020 FP-PSO [34] hybrid / IEEE-69 2020 NLP-PLS MISOCP [40] [41] mathematical mathematical 8 8 IEEE-69 IEEE-69&119 2016 2019 26 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation Real-world problem mathematical formulations are derived under certain assumptions that should reflect and model the problems physical behaviour in as much detail as possible. Even with these assumptions, the solution to large scale power systems is not simple. Methods that propose mathematical programming algorithms are usually unsuitable for large power systems and prone to errors when linearising its non-linear characteristics and uncertainties. It is desirable that a solution to any problem should be a global optimum. However, solutions obtained by mathematical optimisation are not necessarily globally optimal. The reason behind this is usually the modelling complexity and linearisation processes which can be quite challenging. [39] These facts make it difficult for these methods to deal effectively with many power system problems through strict mathematical formulation alone. Furthermore, gradient search, linear programming, dynamic programming, sequential quadratic programming, and nonlinear programming are considered as traditional methods, but none of them provides a solution to complex problem or optimal solution within a reasonable time. [7] It should be mentioned; however, that if the mathematical formulations are reflecting the problems behaviour exactly, even though computation time is rather long, mathematical programming algorithms guarantee a highly superior if not globally optimal result. [41] This paper outlines a new approach to solving the problem of optimal allocation and sizing of PINs by presenting a simple search-based algorithm that in its essence is a pattern of using a load flow procedure. The main motivation behind this research is to introduce a simple solution to a highly complex non-linear problem. By iteratively looping through the network buses, the algorithm places a PIN of a defined size and type (active, reactive, or apparent) at specific locations yielding minimum power losses while keeping all voltages at their acceptable levels. Moreover, the optimal power factor is also determined in the case of apparent type of PIN. The proposed method is demonstrated on a 69-bus and 119-bus distribution system. Table 1 presents an overview of methods from all mentioned groups of optimisation approaches over the past decade, along with several methods/approaches they outperform. Results from the proposed method are compared to those from Tab. 1. It can be concluded that the proposed method is superior to every recent methodology, apart from a very few exceptions. It produces repetitive and unique results; it is extremely easy to implement, and it has a short computation time. 2 PROBLEM FORMULATION Determining the optimal location and size of a PIN in the distribution system is a complex nonlinear problem. As previously mentioned, for a certain location (bus) in the network, as the size of the PIN increases, an adequate decrease in losses can be noticed. However, after exceeding a certain PIN threshold, the losses start to increase again due to a reverse power flow towards the slack/supply bus. The size of the PIN should at most be such that it is consumable within the distribution substation limits. [12] The problem of determining the optimal PIN size and locations is quantified through objective function that can take any mathematical formulation in terms of complexity. For purposes of comparison to other relevant research, objective function introduced in this paper includes costs for energy and power losses only and should yield a least possible value, i.e., minimum: min F= ce ⋅ ∆W + c p ⋅ ∆Pmax (2.1) JET 27 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 with ce being the cost for kilowatt-hour of electricity ($/kWh), ΔW being the electricity losses within the observed period (kWh), cp being the cost per kilowatt for reduction in losses ($/kW) and ΔPmax being the maximum power losses during the same period (kW). To determine the minimum value of the objective function, it is necessary to obtain the minimum value of power losses which would also result in a minimum value of electricity losses during the observed period. The number of constraints and their formulation can be regarded as a way of guiding the optimisation algorithm towards a feasible and possibly optimal solution because not necessarily always these two go hand in hand. Imposing a high number of constraints generally narrows the optimisation algorithm search path and while providing a more accurate and realistic approach, the trade-off is a feasible but not necessarily optimal solution and vice versa. For purposes of comparison to other relevant research, during all calculations, constrains on bus voltage values are imposed by setting an upper and lower bound, [44] resulting in the following constraint: Vmin ≤ Vi ≤ Vmax , for i =1, …, NB (2.2) where Vmin= 0.9 pu, Vmax= 1.1 pu and NB being the number of buses in the given network. 3 LOAD FLOW Aiming for optimal placement and sizing of PIN, it is first and foremost indispensable to determine the parameters of interest in the distribution system, i.e., bus voltages and power losses. Due to distribution systems specific nature such as high R/X ratio and radial structure, conventional load flow methods usually fail to give satisfactory results. [45] Of all proposed power flow solution methods for distribution systems, backward/forward sweep methods are most widely used because of their computational efficiency and robust convergence characteristics. [46-48] The efficiency of the sweeps can be enhanced with oriented branch numbering [49], the only requirement being that the sending bus number i is smaller than the receiving bus number k, i.e., i < k (Fig. 1). Indices from the branch’s sending buses are stored in a vector f, such that i = f(k), where k is the index of the branch’s receiving bus. Additionally, by introducing a fictitious branch with index 1 (sending end index 0), the number of branches NL becomes equal to the number of buses NB, making the sweep procedure very simple and efficient. Figure 1: Branch representation: branch k between buses i (sending end) and k (receiving end) 28 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation Backward sweep consists of equations that calculate the power flow through branches starting from the last one and proceeding in a backward direction towards the supply/slack bus. First, receive receiving end branch flows are calculated using (3.1), where Sdemand is the power demand at the branch’s receiving end and Ykshunt is shunt admittance connected to bus k due to capacitance of the lines and/or connected capacitors (3.2). receive Skreceive = Sdemand + (Ykshunt )∗ ⋅ Vk2 , for k = 1…, NB Ykshunt = j ( Bklines + Bkcap ), for k = 1…, NB (3.1) (3.2) receive receive Sdemand Sdemand − SkPIN , for k = 1…, NB , new = (3.3) Afterwards, sending end branch flows are calculated using (3.4) and added to the branch’s receiving end that supplies them (3.5). send S= k Skreceive + Z kbranch Skreceive Vk 2 ,= for k NL, NL − 1…, 2 (3.4) Sireceive Sireceive + Sksend , = for k NL, NL − 1…, 2 ,= new 3.5) Forward sweep is performed to determine the voltage drops and actual voltages of each bus starting from the slack bus and proceeding in the forward direction towards the last bus using (3.6). Vk = Vi − Z kbranch  S send ⋅ k  Vi  ∗  2, …, NL  , for k =   (3.6) After completing a sweep, the calculated voltages of the present iteration are compared to those from the previous one. If the voltage mismatch between two consecutive iterations is less than the specified tolerance of ε = 10-4, convergence can be achieved. Otherwise, the procedure is repeated until convergence of the solution is attained. After determining the power flow through the branches, it is easy to calculate the active power losses by simply subtracting the real parts from the complex sending and receiving end branch power flows (3.7): Pkreceive (3.7) where Pkreceive is active power from the branch’s sending bus and Pkreceive is active power from the branch’s receiving bus. JET 29 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski 4 JET Volume 16 (2023) p.p. Issue 2, 2023 CLUSTERING OPTIMISATION As has been previously mentioned, many approaches have been proposed in terms of determining the best possible location and size of PINs in an existing distribution system. However, some of them only consider PINs with unity power factor, some are applicable only for single power injecting unit allocation, others offer nearly optimal solutions, and some have inconveniently high computation time. To overcome these shortcomings, a simple search-based clustering optimisation (CO) is presented in this paper. Basic idea behind the CO is by iteratively probing all buses from the distribution system apart from the slack bus, to place PINs of user-defined size and number of locations at locations that would yield least possible losses, which is quantified through (2.1). Power injections can be considered as purely active, purely reactive, or apparent. Iterative bus probing eliminates the need of using any kind of sensitivity analysis for bus selection as the process itself implicitly does that. In case of apparent power injection, the optimal value of the power factor is also determined. Proposed approach does not constrain the power factor value when searching for its optimal due to several reasons: [4-7] • Most if not all utility owners demand from their dispersed generation PINs reactive power support without any specified value or quantity; • Most dispersed generation PINs are owned by private entities which are usually if not exclusively more incentivised when injecting purely active power, i.e., they tend to operate at unity power factor. In some cases, they’re required to operate at unity power factor at their point of common coupling, i.e., they should only produce reactive power for their own personal operational requirements; • Imposing constraints on power factor value implicitly disables the proposed approach in reaching competitive and comparable results to other relevant methods. During the CO procedure, all buses apart from the slack are considered candidate buses for optimal PIN placement, i.e., there is no bus selection procedure nor weak bus sensitivity analysis. There are many reasons for avoiding this type of bus pre-processing. Not all buses possess the same sensitivity trend when subjected to a same rate of power injection change. This is owed to network topology, load, and voltage profile. For example, if a single small power injection is placed on every bus successively, they will attain a certain sensitivity index based on the applied analysis. This index can be used to rank the buses in terms of their susceptibility towards receiving that particular power injection. However, adding a successive power injection of same size introduces a shift in bus sensitivity index following its appropriate trend that is unknown, meaning there may be a shift in position from the previous bus ranking. For a fixed number of locations for PIN placement this is a huge problem since the set of locations potentially changes by successively adding more power injections per bus and the analysis becomes obsolete. There are ways one can alleviate that, i.e., sensitivity analysis and bus ranking can be performed for base case scenario and the obtained set of buses can be kept constant throughout the process of optimisation. However, one should bear in mind that different sensitivity analysis imposes different set of candidate buses and consequently a different solution that may or may not be optimal. [50] 30 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation The proposed approach alleviates this drawback by probing each bus with a small power injection, i.e., implicitly choosing locations that yield least possible losses while inherently following the buses sensitivity trend, hence avoiding the bus selection procedure. The CO algorithm calculates power losses in the distribution network using a linearised (flat start, single iteration) load flow [48] while keeping the voltage profile of the system within prescribed limits. The use of linearised load flow implies that load flow procedure finishes after one iteration, i.e., the resulting power losses and voltages are obtained after a single iteration of backward/forward sweeps. The latter enables placing PINs at each bus successively while keeping the computation time for the entire procedure significantly short. Values of power losses obtained from a linear load flow are smaller compared to those obtained from regular iterative load flow due to flat start, i.e., all bus voltages are equal to the distribution system nominal voltage. Therefore, to obtain the actual values of power losses and voltages, another non-linear load flow is performed at the end of the clustering procedure. This approach is legitimate as the difference in values obtained from both linear and non-linear load flow solution is equal for each bus, therefore the optimal PIN size and their optimal locations do not change. The number of locations for PIN placement, type, size and ϕ step (in case of apparent PIN for attaining optimal power factor) are all user-defined. The CO algorithm performs through the following steps: • Step 1. Read distribution system line and load data and obtain total number of buses Nloc . Initialise user-defined PIN size and type, desired number of locations for PIN placement Nloc , power factor angle step n = 0 and counter n = 0 (current number of buses where PINs are placed). Perform base case load flow and obtain values for distribution system power losses ∆P0 and minimum voltage U min,0 . • ( ) PIN Step 2. Place PIN units with predefined size Sunit and type (active, reactive, or apparent) at each bus ϕ step . Perform subsequent linearised load flows (flat start, single iteration) for each PIN at each bus i, to check for power loss reduction. If apparent PIN type is considered, perform additional linearised load flows for each value of ϕ step to check for optimal power factor as well, for each −π / 2 at each bus. The value of −π / 2 goes from −π / 2 to +π / 2 with a resolution of SiPIN . Stop placing PINs per bus i if no loss reduction and no voltage improvement is achieved. Store the cumulative PIN for bus i as SiPIN and its appropriate losses ∆Pi and minimum voltage Vimin . • Step 3. Using the results from Step 2, for the set of buses= i 2,3, … N B , find bus m that ensures minimum active power losses ∆P= m min( ∆P2 , ∆P3 , …, ∆PN B ) . Place the power injection at SmPIN location m and increase the counter n by 1. JET 31 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski • • JET Volume 16 (2023) p.p. Issue 2, 2023 Step 4. Update the power demand with n = Nloc at the target bus m obtained from Step 3 using (3.3) and repeat the previous two steps until the desired number of locations is achieved ( n = Nloc ) or no power loss reduction is detected. Step 5. Perform final non-linear load flow for the cumulative PINs and placement locations to obtain actual values for power losses and voltages. Figure 2 represents a flow chart of the CO algorithm to further support and explain its performance. Additionally, an example for illustrative purposes is presented below. Example: The CO algorithm is illustrated on a 12.66 kV, 69-bus distribution system searching a single location by placing apparent PINs of size ϕ kVA. The value of ϕ goes from +π / 2 to +π / 2 with ϕ step = π / 9 . The main purpose is power loss minimisation, i.e., minimum value of (2.1) with ∆P0 = 225 kW and lowest voltage 225 . Base case losses are ∆P0 = 225 and ∆P0 = @ 65 occurs at bus 65, i.e., U min,0 = 0.9092 pu. All buses apart from the slack bus are considered for potential PIN placement. Considering the use of apparent PIN in this example, the procedure also determines the optimal value of ϕ step = π / 9 which yields minimum losses, by performing subsequent linearised load PIN flows for each value of ϕ step = π / 9 from +π / 2 to +π / 2 . The latter is done for each Sunit per bus. The CO builds the cluster by continuously adding ϕ per bus and obtaining the optimal value of ϕ until no further loss reduction is detected. When no further power loss reduction is detected by adding additional SiPIN at a certain bus, the cluster, i.e., cumulative PIN with size SiPIN and cumulative value of ϕ is formed. Figure 3 shows a normalised histogram of PINs and power losses for every candidate bus. Losses are normalised to ∆P0 and PINs are normalised to the maximum system injection per bus, which for this example is 4944.7 kVA. 32 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation Figure 2: Flow chart of the CO algorithm JET 33 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 Figure 3: Normalised histograms of PINs and power losses for all candidate buses from the 69-bus distribution system Figure 4: Cluster formation at bus 61 from the 69-bus distribution system PIN Once the cumulative PIN for every bus ( S@ bus ) is determined along with its appropriate losses, the CO performs a simple search for minimum losses and optimal bus for PIN placement. The cumulative PIN is then placed at the determined/optimal location. In this example, the optimal PIN j 35.42 location is bus 61 (Fig. 3), while the cumulative PIN value is = S@ kVA. 61 2173.5 ⋅ e o The cluster formation for the optimal location is visually presented in Fig. 4 and Tab. 2. 34 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation Table 2: Cluster formation at bus 61 from the 69-bus distribution system Iteration № PIN Sunit ϕ 1-6 200ej30 300 7 200ej50 500 8 200ej30 300 9 200ej50 500 10 200ej30 300 11 200ej50 500 2173ej35.42 35.420 Step 4 is omitted in this example since there is only one potential location for PIN placement. @ 61 @ 61 PIN However, the updating of power demand at bus 61 ( Sdemand = , new S demand − S@ 61 ) is still performed considering the new condition in the distribution system. Next, the CO performs a final non-linear load flow to determine the actual power losses @ 27 ∆P = 23.3370 kW and minimum voltage U min = 0.9720 pu in the distribution system. 5 ROBUSTNESS, CONVERGENCE PROPERTIES AND COMPUTATION TIME The proposed CO method possesses several remarkable features, making it superior in terms of complexity/simplicity, implementation, and results compared to other existing search-orientated methods (referenced in Section 1) used for optimal power injection sizing and allocation. • There is no bus selection procedure, nor any form of mathematical modelling and analysis for bus selection. The CO method implicitly chooses locations/buses based on a simple search described in Section 4; • The CO method produces unique, mostly superior, repetitive, and easily reproducible results compared to other search-orientated methods; • The algorithm operates in relatively short computation time that is dependent on distribution system size and PIN type/size. In case of apparent PINs, computation time is generally longer because optimal power factor is also being determined. Moreover, smaller size of PIN units engenders longer computation time due to an increased number of performed calculations per bus; • Users are disburdened from complex mathematical formulations or models as in its essence, the CO method is a ‘pattern of using a load flow procedure’ making it extremely easy to implement. JET 35 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 Table 3 summarises the CO’s performance for the 69-bus distribution system. It can be noted that for rather large range of PIN sizes the power loss and cumulative PIN deviations are small. This is owed to the fact that cluster formation with a certain PIN size can be considered as cluster resolution. Larger PIN size imposes faster and rougher cluster formation, i.e., shorter computation time. Smaller PIN size imposes slower and smoother cluster formation but at the expense of longer computation time. Furthermore, PIN size and type are user-defined variables which can be considered a drawback if one’s so insistent on that tenth or hundredth of a kilowatt losses delta. Real-life distribution system line and load data acquisition can be quite a difficult task. Additionally, distribution system data collection can introduce errors of magnitude in several tens of percent so it can be safely said that the method is as accurate as the accuracy of the input data. Tab. 3 should only serve for user’s illustrative and indicative purposes when choosing PIN size regarding output results and computation time. The CO algorithm is implemented in MATLAB R2018a. All calculations are performed on a laptop configuration with a 2.6-GHz Intel Core i5-4210M processor with 8Gb of RAM. All solutions from references subject to comparison with the CO algorithm are ran through a load flow procedure to check for accuracy of the results. Figure 5 illustrates convergence curves for the CO algorithm when searching for a single location and three locations accordingly. It should be noted that the ending point of the single location convergence curve is a starting point for the next one etc. The latter is owed to the fact that CO algorithm successively chooses buses (one by one) for PIN placement which is explained in detail in Section 4. Table 3: Power loss deviation and calculation time depending on PIN size, type, and number of locations for the 69-bus distribution system Active PIN № PIN Size (kVA) 1 (kW) 36 JET Total (sec.) (kW) ΔPQ (kW) tQ Total ΔPS (sec.) (kVAr) (kW) tS (sec.) Total (kVA) 83.26 8.69 1841 152.06 6.01 1307 23.17 113.13 2237.65ej35.45 50 83.24 0.20 1850 152.08 0.13 1300 23.19 2.23 2223.31ej35.30 100 83.41 0.09 1800 152.08 0.07 1300 23.2 1.15 2273.10ej35.18 200 83.41 0.06 1800 153.43 0.04 1200 23.34 0.60 2173.52ej35.42 71.8 12.24 2353 146.49 8.34 1654 7.57 149.32 2851.33ej35.16 50 71.85 0.25 2350 146.51 0.23 1650 7.51 3.00 2797.55ej34.27 100 71.83 0.19 2300 146.61 0.11 1600 7.88 1.62 2867.80ej34.79 200 72.42 0.10 2300 147.78 0.06 1600 7.5 0.91 2765.32ej35.68 70.28 13.39 3072 145.70 9.99 2168 5.29 159.29 3734.45ej35.24 50 70.34 0.29 3050 145.72 0.25 2150 5.22 3.67 3706.98ej35.03 100 70.32 0.17 3000 145.83 0.17 2100 5.59 1.87 3755.30ej35.23 200 69.78 0.10 2600 146.99 0.10 2200 5.24 1.41 3556.18ej35.52 1 3 tP Apparent PIN 1 1 2 ΔPP Reactive PIN Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation Figure 5: Convergence curves for power losses of the CO method: a) single location 69-bus, b) three locations 69-bus 6 CASE STUDIES Assessment of the effectiveness of the proposed CO method is achieved by its application on two distribution systems: 12.66 kV, 69-bus distribution system [51] and an 11 kV, 119-bus distribution system. [52] To emphasise the unique properties of the CO method, results obtained from the two distribution systems are compared to results from recent work. For purposes of comparison, data concerning the objective function (1) is given with c p = 1 and c p = 1 . Otherwise, objective function values in Tab. 4 and Tab. 6 are calculated with ce = 0.067 $/kWh and c p = 16 $/kW. Base-case values for power losses, minimum voltages, and objective function (1) for both systems are: 69-bus. U min = 0.9092 (kW), U min = 0.9092 (pu) and F0 = 135657 ($/year) 119-bus. ∆P0 = 1298 (kW), F0 = 782590 (pu) and F0 = 782590 ($/year) 6.1 69-bus distribution system Table 4 shows results from the performed analysis using the proposed CO method on the 69bus distribution system. Nine different scenarios are considered, i.e., three scenarios for each PIN type (active – P, reactive – Q and apparent – S) for one, two and three locations. For each scenario, minimum voltages and computation time is also presented. Table 4 suggests that the CO method performs remarkably well presenting unique results in terms of power losses, minimum voltage, and computation time. JET 37 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 Table 4: Results from CO for 69-bus distribution system (power injection size is expressed in P(kW), Q(kVAr) and S(kVA); voltages in pu) PIN № & Type Size@Bus 1P 1869.3@61 2P 1836@61 516@17 3P 1755@61 390@21 1Q 1313@61 2Q 1296@61 336@17 3Q 1302@61 350@17 1S 2249.5ej35.33@61 2S 2209.5e @61 649.6e @17 3S 2120.5e @61 471.4e @21 j35.29 j35.55 Umin@Bus 390@11 ΔP (kW) F ($/ t year) (sec) 0.9683@27 83.22 50175 0.14 0.9807@65 71.77 43272 0.81 0.9795@65 69.7 42024 0.09 0.9305@65 152.05 91674 0.38 j34.00 j35.00 0.9314@65 146.48 88316 0.16 518@50 471.4e j35.00 0.9315@65 145.68 87833 0.57 @11 0.9725@27 23.17 13970 2.22 0.9943@50 7.44 4486 3.64 0.9973@50 4.6 2773 2.90 Table 5: Table of comparison for 69-bus distribution system Method Approach Scenario AM [8] analytical 1P IA [10] analytical AA [13] analytical EA [14] analytical ΔP (kW) t (sec.) CO outperforms 92.00 0.078 Yes 3P 69.96 0.71 Yes 3S 12.72 / Yes 3P 70.20 / Yes 3S 5.92 / Yes 1P 83.23 0.2 Yes 1S 23.26 / Yes 3P 69.70 / Yes 3Q 145.30 / No 152.15 / Yes NH [15] heuristic BPSO-SLFA [26] heuristic 1P TLBO [42] meta-heuristic 3P 72.40 / Yes QOTLBO [42] meta-heuristic 3P 71.62 0.078 Yes IWO [27] meta-heuristic 3P 74.59 5.7 Yes 3S 13.64 / Yes CTLBO [43] meta-heuristic 3P 69.43 / No DE-TLCHS [25] meta-heuristic 3S 5.00 / Yes MRFO [28] meta-heuristic 3P 69.43 / No 38 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation HPSO [29] hybrid MFO [37] hybrid BIBC [36] hybrid QODELFA [32] hybrid LSF-SSA [33] hybrid DGPI [35] hybrid FP-PSO [34] hybrid NLP-PLS [40] mathematical MISOCP [41] mathematical 1P 83.37 / Yes 3S 4.61 / Yes 1P 83.22 / Yes 2S 9.47 / Yes 3P 91.54 / Yes 3P 69.43 / No 3S 12.07 / Yes 1P 83.17 / No 1S 27.91 / Yes 1P 83.22 0.24 Yes 1S 23.30 0.28 Yes 1S 61.67 / Yes 2S 43.67 / Yes 1P 83.23 0.078 Yes 1S 23.18 0.078 Yes 3P 69.43 / No 3S 4.27 / No Table 5 compares the results from the CO method to other methods that use various approaches. Comparison is made for the simplest and the most complicated case presented in the cited references. All results are run through a load flow program to check for accuracy and adequate comparison. CO outperforms in almost every case. It is important to bear in mind that power loss differences per scenario, amongst all compared methods, differ within tenth of a kilowatt at best and several kilowatts at worst, regardless the method and its approach. Given the uncertainty of real distribution systems input data, it can be said that methods accuracy strongly depends on the accuracy of the input data. This is where CO strongly outperforms any other method as it is extremely simple for implementation - pattern of using a load flow procedure, something that every power system engineer knows extremely well. Not everyone possesses a strong modelling and mathematical knowledge and background to shape these types of problems, but certainly everyone or almost everyone knows how to implement a load flow procedure and do a simple search. In terms of computation time (wherever noted, ‘/’ in Table 5 means unreported), CO outperforms most of the methods. Mathematical methods depending on the problem size could take several tens of minutes or even hours in some cases. [41] Analytic, heuristic and hybrid methods possess computation time comparable to the CO method. 6.2 119-bus distribution system Performance evaluation of the CO method for the 119-bus distribution system is attained through five scenarios including: three scenarios for optimal placement and sizing of active, reactive, and apparent PINs at five locations, and two more scenarios including optimal placement and sizing of active and apparent PINs with optimal power factor at seven locations. JET 39 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 Table 6 shows results from the performed analysis using the CO method on the 119-bus distribution system. Comparison results are presented in Tab. 7. Again, all results are run through a load flow program to check for accuracy and adequate comparison. CO outperforms in almost every case, except when compared to MISOCP. The reason behind this is that in [41], an exact model for the distribution system and load flow equations is used with no linearisation and no neglections, meaning that the obtained results present a global optimum. This is a mathematical approach and takes quite a while to reach an optimal solution (several minutes to several hours depending on the scenario). However, the CO method gives results that are very near those presented with MISOCP and in reasonable computation time as it can be seen from Tab. 6. Compared to other methods, CO outperforms in terms of power losses, i.e., solutions differ from several kilowatts at best to several tens of kilowatts at worst in favour of the proposed algorithm. Computation time is better, if not comparable, to other methods subject to comparison. Table 6: Results from CO for 119-bus distribution system (power injection size is expressed in P(kW), Q(kVAr) and S(kVA); voltages in pu) Size@Bus PIN № & Type Total PIN (kVA) ΔP (kW) Umin (pu) F ($/year) t (sec) 5P 5Q 5S 7S 7P 2800@50 2800@71 2400@79 1600@96 2800@110 4400@29 2600@50 1800@72 1700@80 2300@110 3942ej42.0@50 3179ej32.5@72 2765ej35.7@80 1981ej34.0@96 3546ej38.9@110 1760@20 1870@41 2860@50 2860@71 2200@80 1540@97 3080@109 2072ej36.6@20 2420ej35.7@41 3796ej42.8@50 3131ej32.2@72 2760ej37.5@80 2082ej33.3@96 3447ej40.0@110 12400 574.33 0.955 346280 0.82 12800 857.69 0.9074 517120 2.14 15413 211.58 0.9608 127570 9.42 16170 515.82 0.9582 311000 3.13 19708 127.61 0.9764 76940 13.08 Table 7: Table of comparison for 119-bus distribution system Method Approach NH [15] heuristic QOTLBO [42] meta-heuristic 40 JET Scenario ΔP (kW) t (sec) CO outperforms 5P 5Q 5S 7P 7S 7P 580.74 861.60 227.50 515.70 128.80 576.00 2.5 2.0 5.4 4.1 6.2 20.83 Yes Yes Yes No Yes Yes Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation CTLBO [43] meta-heuristic LSF-BFOA [31] hybrid QODELFA [32] hybrid MISOCP [41] mathematical 7 7P 5P 5Q 5S 7P 7S 7P 7S 5P 5Q 5S 7P 7S 516.25 578.97 871.40 227.90 526.34 132.10 518.65 132.79 571.29 856.37 208.13 513.27 123.37 / 23.24 23.21 24.65 24.96 25.96 / / / / / / / Yes Yes Yes Yes Yes Yes Yes Yes No No No No No CONCLUSION This paper presents a clustering optimisation (CO) method for optimal placement and sizing of DG’s presented as active, reactive, or apparent power injections (PINs) with the aim of obtaining minimum power losses while maintaining an acceptable system voltage profile. The effectiveness of the proposed method is assessed through simulations performed on a 69-bus and 119-bus distribution systems. It can be concluded that the proposed method performs outstandingly, presenting unique results in all considered scenarios. Results attained from the performed analysis are compared to results from recent methodologies stretching back over the past decade. Tables 5 and 7, coupled with the information from Tab. 1, suggest that the proposed method is implicitly compared to over 50 other methods. It is noted that for most scenarios, CO gives lower power losses and total PIN size compared to other methods, while for other scenarios it presents slightly higher power losses and total PIN size. Power losses and total PIN sizes attained from the CO method are of the same order of magnitude as power losses and total PIN sizes attained from recent methodologies in all considered cases. Therefore, the slight deviation (fractions of kilowatts) in values obtained from the CO method compared to other methodologies is practically insignificant/negligible, considering the uncertainty of input data. In its essence, the CO method is a pattern of using a load flow procedure. It does not include solving complex mathematical formulations, and it does not operate with population, nor its creation and iterative updating. In addition to this, it does not use problem coding and/or solution decoding and has no convergence problems, the latter of which renders the method superior in terms of its simplicity/complexity and easy implementation. Furthermore, the lack of laborious procedures makes the method significantly intelligible, simple, and easy to reproduce. Since the methods working principle is based on simple mathematical formulations and there is no dedicated bus selection procedure, optimal or nearly optimal solutions are obtained remarkably quickly. Short computation time enhances the convenience of using this method. Finally, another salient advantage of the proposed method is the recurrence of results, i.e., the results obtained are repetitive from simulation to simulation, which is not the case for some of the other existing methods. JET 41 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski 8 JET Volume 16 (2023) p.p. Issue 2, 2023 References [1] T. Ackermann, G. Andersson, L. Soder: Distributed generation: a definition, Electric Power Systems Research 57 (3) (2001) 195 – 204. doi: https://doi.org/10.1016/S03787796(01)00101-8 [2] T. Griffin, K. Tomsovic, D. Secrest, A. Law: Placement of dispersed generation systems for reduced losses, in: Proceedings of the 33rd Annual Hawaii International Conference on System Sciences, 2000, pp. 9 pp.–. doi:10.1109/HICSS.2000.926773 [3] N. S. Rau, Y. H. Wan: Optimum location of resources in distributed planning, IEEE Transactions on Power Systems 9 (4) (1994) 2014–2020. doi:10.1109/59.331463 [4] M. Aman, G. Jasmon, A. Bakar, H. Mokhlis, M. Karimi: Optimum shunt capacitor placement in distribution system - a review and comparative study, Renewable and Sustainable Energy Reviews 30 (2014) 429 – 439. doi: https://doi.org/10.1016/j.rser.2013.10.002 [5] M. Pesaran, P. D. Huy, V. K. Ramachandaramurthy: A review of the optimal allocation of distributed generation: Objectives, constraints, methods, and algorithms, Renewable and Sustainable Energy Reviews 75 (2017) 293 – 312. doi: https://doi.org/10.1016/j. rser.2016.10.071 [6] A. Ehsan, Q. Yang: Optimal integration and planning of renewable distributed generation in the power distribution networks: A review of analytical techniques, Applied Energy 210 (2018) 44 – 59. doi: https://doi.org/10.1016/j.apenergy.2017.10.106 [7] Z. Abdmouleh, A. Gastli, L. Ben-Brahim, M. Haouari, N. Al-Emadi: Review of optimization techniques applied for the integration of distributed generation from renewable energy sources, Renewable Energy 113 (05 2017). doi: 10.1016/j.renene.2017.05.087 [8] T. Gozel, M. H. Hocaoglu: An analytical method for the sizing and siting of distributed generators in radial systems, Electric Power Systems Research 79 (6) (2009) 912 – 918. doi: https://doi. org/10.1016/j.epsr.2008.12.007 [9] D. Q. Hung, N. Mithulananthan, R. Bansal: An optimal investment planning framework for multiple distributed generation units in industrial distribution systems, Applied Energy 124 (2014) 62 – 72. doi: https://doi.org/10.1016/j.apenergy.2014.03.005 [10] D. Q. Hung, N. Mithulananthan: Multiple distributed generator placement in primary distribution networks for loss reduction, IEEE Transactions on Industrial Electronics 60 (4) (2013) 1700–1708. doi:10.1109/TIE.2011.2112316 [11] H. Khan, M. A. Choudhry: Implementation of distributed generation (IDG) algorithm for performance enhancement of distribution feeder under extreme load growth, International Journal of Electrical Power & Energy Systems 32 (9) (2010) 985 – 997. doi: https://doi. org/10.1016/j.ijepes.2010.02.006 [12] N. Acharya, P. Mahat, N. Mithulananthan: An analytical approach for DG allocation in primary distribution network, International Journal of Electrical Power & Energy Systems 28 (10) (2006) 669 – 678. doi: https://doi.org/10.1016/j.ijepes.2006.02.013 42 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation [13] S. N. G. Naik, D. K. Khatod, M. P. Sharma: Analytical approach for optimal siting and sizing of distributed generation in radial distribution networks, IET Generation, Transmission & Distribution 9 (2015) 209–220(11) [14] K. Mahmoud, N. Yorino, A. Ahmed: Optimal distributed generation allocation in distribution systems for loss minimization, IEEE Transactions on Power Systems 31 (2) (2016) 960–969. doi: 10.1109/TPWRS.2015.2418333 [15] A. Bayat, A. Bagheri: Optimal active and reactive power allocation in distribution networks using a novel heuristic approach, Applied Energy 233-234 (2019) 71 – 85. doi: https://doi. org/10.1016/j.apenergy.2018.10.030 [16] K. Muthukumar, S. Jayalalitha: Optimal placement and sizing of distributed generators and shunt capacitors for power loss minimization in radial distribution networks using hybrid heuristic search optimization technique, International Journal of Electrical Power & Energy Systems 78 (2016) 299 – 319. doi: https://doi.org/10.1016/j.ijepes.2015.11.019 [17] M. H. Moradi, A. Zeinalzadeh, Y. Mohammadi, M. Abedini: An efficient hybrid method for solving the optimal sitting and sizing problem of DG and shunt capacitor banks simultaneously based on imperialist competitive algorithm and genetic algorithm, International Journal of Electrical Power & Energy Systems 54 (2014) 101 – 111. doi: https://doi.org/10.1016/j. ijepes.2013.06.023 [18] A. Bayat, A. Bagheri, R. Noroozian: Optimal siting and sizing of distributed generation accompanied by reconfiguration of distribution networks for maximum loss reduction by using a new UVDA-based heuristic method, International Journal of Electrical Power & Energy Systems 77 (2016) 360 – 371. doi: https://doi.org/10.1016/j.ijepes.2015.11.039 [19] Y. Li, B. Feng, G. Li, J. Qi, D. Zhao, Y. Mu: Optimal distributed generation planning in active distribution networks considering integration of energy storage, Applied Energy 210 (2018) 1073 – 1081. doi: https://doi.org/10.1016/j.apenergy.2017.08.008 [20] M. Dixit, P. Kundu, H. R. Jariwala: Incorporation of distributed generation and shunt capacitor in radial distribution system for techno-economic benefits, Engineering Science and Technology, an International Journal 20 (2) (2017) 482 – 493. doi: https://doi. org/10.1016/j.jestch.2017.01.003 [21] N. Kanwar, N. Gupta, K. Niazi, A. Swarnkar, R. Bansal: Simultaneous allocation of distributed energy resource using improved particle swarm optimization, Applied Energy 185 (2017) 1684 – 1693, clean, Efficient and Affordable Energy for a Sustainable Future. doi: https://doi.org/10.1016/j.apenergy.2016.01.093 [22] U. Sultana, A. B. Khairuddin, A. Mokhtar, N. Zareen, B. Sultana: Grey wolf optimizer-based placement and sizing of multiple distributed generation in the distribution system, Energy 111(2016) 525 – 536. doi: https://doi.org/10.1016/j.energy.2016.05.128 [23] A. Zeinalzadeh, Y. Mohammadi, M. Moradi: Optimal multi objective placement and sizing of multiple DGs and shunt capacitor banks simultaneously considering load uncertainty via mopso approach, International Journal of Electrical Power & Energy Systems 67 (05 2015). doi: 10.1016/j.ijepes.2014.12.010 JET 43 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 [24] M. Moradi, M. Abedini: A combination of genetic algorithm and particle swarm optimization for optimal DG location and sizing in distribution systems, International Journal of Electrical Power & Energy Systems 34 (1) (2012) 66 – 74. doi: https://doi.org/10.1016/j. ijepes.2011.08.023 [25] M. Zaid, A. Alam, Z. Sarwer, U. Shahjahani: Optimal allocation of distributed energy resources in a distribution system, in: Conference: INTERNATIONAL CONFERENCE ON COMPUTING, POWER AND COMMUNICATION TECHNOLOGIES 2018 (GUCON) At: Galgotias University, Noida, 2019. doi:10.1109/i-PACT44901.2019.8959995 [26] A. Shuaibu, Y. Sun, Z. Wang: Multi-objective for optimal placement and sizing DG units in reducing loss of power and enhancing voltage profile using BPSO-SLFA, Energy Reports 6 (2020) 1581–1589. doi: 10.1016/j.egyr.2020.06.013 [27] D. Prabha, T. Jayabarathi: Optimal placement and sizing of multiple distributed generating units in distribution networks by invasive weed optimization algorithm, Ain Shams Engineering Journal 7 (08 2015). doi: 10.1016/j.asej.2015.05.014 [28] M. G. Hemeida, A. A. Ibrahim, A.-A. A. Mohamed, S. Alkhalaf, A. M. B. El-Dine: Optimal allocation of distributed generators DG based manta ray foraging optimization algorithm (MRFO), Ain Shams Engineering Journal (2020). doi: https://doi.org/10.1016/j. asej.2020.07.009 [29] S. Kansal, V. Kumar, B. Tyagi: Hybrid approach for optimal placement of multiple DGs of multiple types in distribution networks, International Journal of Electrical Power & Energy Systems 75(2016) 226 – 235. doi: https://doi.org/10.1016/j.ijepes.2015.09.002 [30] M. Gitizadeh, A. A. Vahed, J. Aghaei: Multistage distribution system expansion planning considering distributed generation using hybrid evolutionary algorithms, Applied Energy 101 (2013) 655 – 666, sustainable Development of Energy, Water and Environment Systems. doi: https://doi.org/10.1016/j.apenergy.2012.07.010 [31] K. Devabalaji, K. Ravi: Optimal size, and siting of multiple DG and DSTATCOM in radial distribution system using bacterial foraging optimization algorithm, Ain Shams Engineering Journal 7 (3) (2016) 959 – 971. doi: https://doi.org/10.1016/j.asej.2015.07.002 [32] R. Mahfoud, Y. Sun, N. Alkayem, H. Haes Alhelou, P. Siano, M. Shafie-khah: A novel combined evolutionary algorithm for optimal planning of distributed generators in radial distribution systems, Applied Sciences 9 (08 2019). doi:10.3390/app9163394 [33] M. J. Tahir, B. A. Bakar, M. Alam, M. M. Su’ud: An LSF based DG placement approach for power losses minimization of the radial distribution network, in 2019 13th International Conference on Mathematics, Actuarial Science, Computer Science and Statistics (MACS), 2019, pp. 1–6. doi:10.1109/MACS48846.2019.9024811 [34] S. Kumar, S. K. Ra: Optimal integration of distributed generators using hybrid flower pollination algorithm, IJEET (08 2020). doi:10.13140/RG.2.2.14868.12169 [35] G. Memarzadeh, F. Keynia: A new index-based method for optimal DG placement in distribution networks, Engineering Reports 2 (07 2020). doi:10.1002/eng2.12243 44 JET Optimal Dispersed Generation Placement in Radial Distribution Systems Using a Pattern of Load Flow Procedure – Clustering Optimisation [36] A. Alam, M. Zaid, A. Gupta, P. Bindal, A. Siddiqui: Power loss reduction in a radial distribution network using distributed generation, in 2018 International Conference on Computing, Power and Communication Technologies (GUCON), 2018, pp. 1142–1145. doi:10.1109/GUCON.2018.8674942 [37] H. Abdel-Mawgoud, S. Kamel, M. Ebeed, A. Youssef: Optimal allocation of renewable DG sources in distribution networks considering load growth, in 2017 Nineteenth International Middle East Power Systems Conference (MEPCON), 2017, pp. 1236–1241. doi:10.1109/ MEPCON.2017.8301340 [38] S. Nojavan, M. Jalali, K. Zare: Optimal allocation of capacitors in radial/mesh distribution systems using mixed integer nonlinear programming approach, Electric Power Systems Research 107 (2014) 119 – 124. doi: https://doi.org/10.1016/j.epsr.2013.09.019 [39] A. Rueda-Medina, J. Franco, M. J. Rider, A. Padilha-Feltrin, R. Romero: A mixed-integer linear programming approach for optimal type, size, and allocation of distributed generation in radial distribution systems, Electric Power Systems Research 97 (2013) 133– 143. doi: 10.1016/j.epsr.2012.12.009 [40] V. Pamishetti, S. Singh: Optimal placement of dg in distribution network for power loss minimization using NLP & PLS technique, Energy Procedia 90 (2016) 441–454. doi: 10.1016/j.egypro.2016.11.211 [41] J. Vuletic, J. Angelov and M. Todorovski: Optimal power injection placement in radial distribution systems using mixed integer second order cone programming, in: International Scientific Conference on Information, Communication and Energy Systems and Technologies (ICEST) 2019, june, 2019 [42] S. Sultana, P. K. Roy: Multi-objective quasi-oppositional teaching learning-based optimization for optimal location of distributed generator in radial distribution systems, International Journal of Electrical Power & Energy Systems 63 (2014) 534 – 545. doi: https://doi.org/10.1016/j.ijepes.2014.06.031 [43] A. Quadri, S. Bhowmick, D. Joshi: A comprehensive technique for optimal allocation of distributed energy resources in radial distribution systems, Applied Energy 211 (2018) 1245 – 1260. doi: https://doi.org/10.1016/j.apenergy.2017.11.108 [44] EN 50160 - voltage characteristics of electricity supplied by public electricity networks. [45] S. C. Tripathy, G. D. Prasad, O. P. Malik, G. S. Hope: Load-flow solutions for ill-conditioned power systems by a newton-like method, IEEE Transactions on Power Apparatus and Systems PAS-101 (10) (1982) 3648–3657. doi:10.1109/TPAS.1982.317050 [46] D. Shirmohammadi, H. W. Hong, A. Semlyen, G. X. Luo: A compensation-based power flow method for weakly meshed distribution and transmission networks, IEEE Transactions on Power Systems 3 (2) (1988) 753–762. doi:10.1109/59.192932 [47] G. X. Luo, A. Semlyen: Efficient load flow for large weakly meshed networks, IEEE Transactions on Power Systems 5 (4) (1990) 1309–1316. doi:10.1109/59.99382 JET 45 Angela Popova, Jovica Vuletić, Jordančo Angelov, Mirko Todorovski JET Volume 16 (2023) p.p. Issue 2, 2023 [48] M. E. Baran, F. F.Wu: Network reconfiguration in distribution systems for loss reduction and load balancing, IEEE Transactions on Power Delivery 4 (2) (1989) 1401–1407. doi:10.1109/61.25627 [49] D. Rajicic, R. Ackovski, R. Taleski: Voltage correction power flow, IEEE Transactions on Power Delivery 9 (2) (1994) 1056–1062. doi:10.1109/61.296308 [50] S. Raj, B. Bhattacharyya: Weak bus-oriented optimal var planning based on grey wolf optimization, in 2016 National Power Systems Conference (NPSC), 2016, pp. 1–5. doi:10.1109/NPSC.2016.7858929 [51] D. Das: Optimal placement of capacitors in radial distribution system using a Fuzzy-GA method, Int. J. Elect. Power Energy Syst. 30 (6-7) (2008) 361 – 367. doi: http://dx.doi. org/10.1016/j.ijepes.2007.08.004 [52] S. Ghasemi: Balanced and unbalanced distribution networks reconfiguration considering reliability indices, Ain Shams Engineering Journal 9 (4) (2018) 1567 – 1579. doi: https://doi. org/10.1016/j.asej.2016.11.010 46 JET JET Volume 16 (2023) p.p. 47-56 Issue 2, 2023 Type of article: 1.01 http://www.fe.um.si/si/jet.htm VISUALISATION OF 4D THERMAL MAPS 4D VIZUALIZACIJA TEMPERATURNEGA POLJA Gorazd Hren R1 Keywords: real-time visualisation, volume rendering, spatial interpolation Abstract Today, there are many building energy simulation systems. Most rely on theoretical thermal models and numerical simulations to improve the building design and achieve good energyefficient living conditions. The 3D visualisation technique proposed in this paper analyses temperature distribution in a room. A room is divided into uniform grids, and the temperature field is established by calculating temperature values for each cell using an inverse distance weighted interpolation algorithm with sensor data. These cells are visualised by voxel rendering. An X3D simulation system coupled with a cost-effective set of temperature sensors can provide valuable insights into room design. A measured thermal status can lead to energy savings and improved thermal comfort with improved efficiency. The paper discusses the cost-effective and simplified process for a hardware-software prototype system that can provide real data-driven 4D thermal maps for buildings. The system can be extended to provide 4D thermal maps for more rooms and multi-storey buildings. Povzetek Dandanes obstaja veliko aplikacij za energijsko simulacijo stavb in večina uporablja teoretične toplotne modele in numerične simulacije za izboljšanje zasnove stavbe in doseganje energetsko učinkovitih pogojev človeškega udobja. V tem prispevku je za analizo porazdelitve temperature v prostoru predlagana tehnika 3D vizualizacije v predpisanih časovnih okvirjih. Prostor je razdeljen na mrežo celic, temperaturno polje je izračunano z uporabo algoritma interpolacije inverzne razdalje. Iz podatkov o poziciji in vrednosti temperature v posameznem R Corresponding author: Assoc. Prof. Gorazd Hren Ph. D., University of Maribor, Faculty of Energy Technology, Hočevarjev trg 1, 8270 Krško, Gorazd.Hren@um.si JET 47 Gorazd Hren JET Volume 16 (2023) p.p. Issue 2, 2023 senzorju v prostoru interpoliramo vrednosti temperaturnega polja v celotnem prostoru. Za vizualizacijo uporabimo metodo upodabljanja volumetričnih slikovnih pik (voxel). Simulacijski sistem X3D skupaj s stroškovno učinkovitim naborom temperaturnih senzorjev lahko zagotovi vpogled v temperaturno polje prostora na osnovi izmerjenih temperatur. Temperaturno polje prikaže stanje temperature v celotnem prostoru v časovnih korakih, kar lahko privede do prihrankov energije, izboljšanega toplotnega udobja z izboljšano učinkovitost upravljanja klimatizacije prostora. Prispevek predstavlja učinkovit postopek za izdelavo 4D temperaturnega polja na podlagi izmerjenih podatkov. Rešitev je mogoče razširiti za časovno vizualizacijo temperaturnega polja za večnadstropne zgradbe. 1 INTRODUCTION Thermal management is an effective method to reduce energy consumption. [1] Thermal management protocols can optimise temperature distribution, decrease heating and cooling costs and increase human comfort. The crucial task of thermal management is understanding the airflow and finding the hot or cold spots. Visualisation of temperature distribution as a 3D temperature field is an efficient method to achieve this. Many attempts only visualise the temperature field in 2D, which then needs to be revised to present the entire temperature in space. [2,3] For the most part, computational fluid dynamics (CFD) methods, [4] which are time consuming and neglect the actual temperature state in spaces, lake home rooms, warehouses, or production halls, are adopted. However, its effectiveness is limited by the computational costs associated with predictions, which limit online operation and maintenance applications. Consequently, effective temperature assessments still depend on the limited measurement data obtained by discrete temperature sensors. [5,6] Nowadays, many closed spaces are grided with temperature sensors. A method to produce 3D visualisation for thermal management with the temperature distribution by 3D particles from temperature data in time extracted from sensors, [7,8] but the results still need to be understood. This approach first establishes the temperature field with data from sensors and, after visualising the temperature distribution field in time, is performed. Three aspects are taken into consideration to gain insight into temperature distribution in closed spaces: • Practical 3D visualisation program by modelling temperature as a data volume set. From this, temperature distribution could be visualised as volume rendering. • For calculating the temperature field, the inverse distance weighted interpolation (IDW) is used. Heat transfer rules and transfer distance are used for the method, which makes the visualisation results more reliable. The execution of complex algorithms is a question of the available CPU RAM. [9] • The design rendering method (and variations) for different analysis purposes. Rendering temperature field grid of voxels where each is represented as a box using transparent colour as a visual element. • The approach can be roughly divided into two steps: establishing the temperature field from sensor data and visualising the temperature field. 48 JET Visualization of 4d thermal maps 2 SPATIAL TEMPERATURE FIELD CALCULATION The first step is to divide the room into a uniform grid, with each cell becoming the element for the visual component. The temperature of each cell is calculated by the IDW method, which simulates the sensors as heat sources and bounds the transfer distance of sensors. IDW reduces the complexity of computation and maintains the locality. The temperature field is divided into space voxels, transforming the problem into volume data. [10-13] The abstracted model is a single room r=(w,h,d), w-width, h-height, and d-depth. The room is divided into a 3D grid that is regularly spaced. The cell size determines the density of the grid and space division. Each cell is considered a voxel, and the whole room represents a collection of voxels. The temperature field describes the temperature distribution in the entire room. The voxels define it, and each cell’s temperature value is defined. The collection of all voxels with temperature data is volume data. The temperature field is equivalent to a set of volume data, and visualising the temperature distribution can be transformed into a volume rendering problem, specifically to calculate colours representing each voxel in space. For visualising the volume data set, the temperature value of each voxel must first be calculated by a spatial interpolation method with data from sensors. Temperature distribution is a result of heat transfer. Voxels close to each other have similar temperature values, while voxels at larger distances are likely different. Thermal distribution is thus distance-related. Inverse distance interpolation IDW is a widely used method with known values in scattered points corresponding to the thermal distribution characteristic. The basic IDW function is shown below: n (2.1) where T is the interpolation value at point x from data Ti=T(xi) and i {0,1,…n}; xi is the known point, d is the Euclidean distance between the estimated point and the observation point, x, and xi, and s is the number of points with known values, and n is the total number of voxels. p is the positive real number, a power parameter; greater values of p assign greater influence on values closest to the interpolation point. The implicit assumption of the classical IDW interpolation method is that the objects to be estimated are uniformly distributed in space with the power parameter p=2. 2.1 Modelling and visualisation in application The abstracted model, a single room 3x3 m (SR33), was produced as a testbed to develop an approach, check different options, and verify the results. MatLab was used as the programming language for calculations and preliminary visualisation. MatLab applies interp1, interp2, and interp3 functions in [14] for uniform data. These functions use polynomial techniques, fitting the supplied data with polynomial functions between data points and evaluating the appropriate function at the desired interpolation points. When the sample data is scattered, the interpolation techniques use a triangulation-based approach to compute interpolated values. The interpolating function passes through the data points. This is an important distinction between interpolation JET 49 Gorazd Hren JET Volume 16 (2023) p.p. Issue 2, 2023 and curve/surface fitting. The function does not necessarily pass through the sample data points in the fitting. We do have measurements in space: the exact position and temperature value at that point. The SR33 is grided with evenly spaced points. The sensors’ location and temperature data are prescribed. In the SR33 corner, a temperature of -10 C is prescribed. In the opposite corner, a temperature of 10 C is prescribed, while in the middle of the SR33, 25 C is prescribed. The temperature field of SR33 was organised as a set of volume data, thus allowing the temperature field to be visualised by rendering methods. Each voxel is assumed to be an element (sphere or box) and is rendered as an element. The colour and scale of transparent elements are changed to determine heat distribution in the room. Visualisation of spatial data is challenging to achieve. In order to visualise the volumetric data set, the IWD must be interpolated from the sensor data to calculate the temperature value of each voxel. The space is modelled with a rough grid for easier and faster calculations. The sensors’ position, their temperature data, and the calculated temperatures in all interpolation points are determined. For the refined mesh, the mapping is performed using the MatLab interp3 spatial interpolation function. The representation and visualisation of room temperature in MatLab is limited, although MatLab contains many functions in the libraries for drawing graphs and pictures. Various functions can be applied to 3D volume data; Figures 1-3 show the various functions and some of their variations and options. Only three functions were relevant: scatter, slice, [14] and VOXview. [15] Figure 1 shows the results using the scatter and VOXview functions, while Figure 2 uses the slice function, which produces multiple slices. Figure 3 presents different IWD interpolations with power parameters p=1 and p=2. In all images, the temperature mapping is scaled between the highest and lowest temperature from all temperature data values. Figure 1: Results of IDW interpolation presented by scatter and VOXview function; with basic grid and with inerp3 refined grid, respectively Visualisation of the temperature field helps to analyse and manage heat flows. Relevant questions include whether there is a hot spot and the overall thermal condition. An overall view and details are required to determine the general thermal state of the room. To produce a visualisation of changes in the temperature field over time, looping with a pause is necessary. There needs to be more quality visualisation since the figures present the temperature field as a sequence over time that cannot be controlled. A detailed view could be achieved with crosssection views, but phenomena such as hot spots or sudden temperature changes could not be predicted in advance. 50 JET Visualization of 4d thermal maps Figure 2: Results of IDW interpolation presented by cross sections and slice function; with basic grid and with inerp3 refined grid, respectively Figure 3: Results of IDW interpolation presented by cross sections and slice function with power parameters of p=1 and p=2; with basic grid and with inerp3 refined grid, respectively 2.2 X3D visualisation X3D is an open standard with worldwide validity certified by the International Organisation for Standardisation. The openness of X3D ensures the independence of a data format. So, like an XML document, an X3D document consists of nodes organised into parent and child elements. Each element has a unique name, defined in the XML Schema Definition. These names are selfexplanatory, and are described in a computer language that is also readable and writeable by users. Four predefined positions are accessible from the view bar at the bottom of the browser: top, front, back, and front-top view. Navigation is defined as both flying and walking through the scene. X3D offers a wide range of functionalities and node types to include data from many sources. There are various ways to create 3D objects in X3D. X3D geometric primitives (e.g., Box, Sphere) are good possibilities for 3D thermal map representation as their RGB and g-values can be controlled straightforwardly, the author chooses to model the volume in the room as a set (mesh) of semi-transparent colour cubes. JET 51 Gorazd Hren JET Volume 16 (2023) p.p. Issue 2, 2023 The cube’s colour represents the temperature at a certain location. Values are interpolated to determine intermediary temperature values. IDW interpolation is used currently, as illustrated in the previous chapter. However, other more complex models may enhance the visual representation of the sensor data. X3D scene graph structure consists of nodes with properties called fields. A field includes an attribute storing the basic data type (e.g. a floating-point value or a vector of three values). The routing of an output field of a geometrical child node of an X3D to the input field of another makes animation within the scene graph possible. X3D event animation at author-time is an event cascade across sensor, interpolator, and transform nodes. The ROUTE mechanism connects the node sending the event and the node receiving that event. A TouchSensor is placed on the keyword geometry used as a trigger. The event is first routed to a TimeSensor node, where the cycle interval is controlled and can be easily modified. To complete the colouring effect, the startTime field of this TimeSensor is routed to a ColorInterpolator, and then to the Material node. 3D worlds can be built without programming, just by describing the scene graph using X3D syntax, particularly in geoscience, BIM. [19] In our case, geometry is not an issue, but it is necessary to transform a vast amount of data, which makes it possible to perform visualisation with computer programming. There has been much experience with visualisations using X3D in the scientific community. A voxel representation algorithm is used to describe the 3D data. This algorithm describes 3D objects with flexible volumes for the whole internal space with predefined shapes (box and sphere) as a 3D mesh. The result of the study should give a good indication of the possibility of producing sufficient 4D maps that could be published via Internet Standards. While X3D is an open standard, powerful, extensible, and has a consortium that maintains and supports it, it needs to provide a more complex environment of quick development. [16] 3 RESULTS AND DISCUSSION The proposed method has been deployed in one room of the Faculty of Energy Technology (LABB) laboratory building. The LABB is a standalone building consisting only of one room with computer hardware equipment. Data collected from 24 sensors are spread in all corners of the LABB (inside, in-wall, and outside) and the middle of the roof. The distribution and temperature values of these sensors are demonstrated in Figure 4. All sensors are linked to a SCADA system that continuously stores data every 5 minutes. The positions of the sensors can be easily seen in Figure 4. The data were obtained in an MS Excel file and reduced to include the temperature values for each hour of the day for each month. Only in-room temperature measurements were considered for visualisation. Since the MatLab application was developed to read and calculate the temperature distribution field in the room, it was possible to change the time step and duration of the measurements by editing the input file. Given the number of sensors deployed in the room is s, and the sensor set is represented by S, then a sensor is denoted by Si, iÎ{1,2, …s} and Si=(x,y,z,T), where values represent the location coordinate in width, height, depth, and the value of the temperature respectively. 52 JET Visualization of 4d thermal maps Figure 4: SCADA layout showing temperature measuring positions sensors The starting point is the LABB geometric model, modelled in the SolidWorks CAD system, extracted to VRML97, and then translated to the currently active X3D format. The walls and glass in the windows are set to a semi-transparent colour, and the computer and projector parts are pink, as can be seen in Figure 5, for easy monitoring. The geometry model and grid of points from the application are linked (the coordinate systems were very helpful). The application reads the temperature values from the Excel file for nine sensors. All considered sensors have a temperature value for each hour for the whole month (720 values for each sensor). Adding more data could significantly slow down the calculation speed. The application transfers data to nodes and X3D files, simultaneously produces all geometric data and all ColorInterpolator data for time frames and calculates a representative colour for the current temperature in a fixed colour bar that ranges from the lowest to the highest temperature value. The results are difficult to show on screen or paper. The temperature in the room is very consistent, with only very small differences between the data of the individual sensors being recorded. The temperature management of the LABB is very efficient. Therefore, two sequences are chosen to show the results: in Figure 6, when the central computer is active, which also has a sensor, and in Figure 7, which shows the temperature field when the projector is active in global view and local detail. JET 53 Gorazd Hren JET Volume 16 (2023) p.p. Issue 2, 2023 Figure 5: CAD model of the LABB and its presentation in X3D viewer (the coordinate system is added for checking data) Figure 6: Temperature distribution field in LABB, on X3D viewer while the computer runs Figure 7: Temperature distribution field in LABB, on X3D viewer, while the projector is on 54 JET Visualization of 4d thermal maps 4 CONCLUSIONS AND FUTURE WORK This paper has presented a basic and cost-effective system for temperature data acquisition and an X3D module for 3D Thermal Map representation. It models the temperature field as a volume data set. The voxel value is obtained by interpolating the sensor data using the IDW method. Voxel rendering is used for visualisation so that the results give a more comprehensive understanding of heat distribution. The proposed method was applied in an actual work environment. The results proved that our method could visualise heat fluxes in 3D and in real time, and help in efficient heat distribution analysis. Since voxel computation and rendering could be done in parallel, the method can be easily implemented on GPU to speed it up. This method can also be used in other applications, such as in the visualisation of humidity or in other environmental factors in closed spaces. In addition to this, optimised sensor distribution will provide better interpolation results. Even though the intelligence is lost due to the transfer of these files (mainly parametric objects) into the X3D environment, the proposed system replaces the virtual model of a building with a three-dimensional virtual prototype containing a swarm of sensors that generates important thermal information ready to be used. As a result, embedded information can be conveyed and input characteristics exported to other analysis tools, such as CFD or ultrasonic method. [17] The most exciting part of this research is the low cost of the system to implement and the visualisation and production of heat maps without specialised software. The X3D concept of paths, interpolators, and sensors for animating elements still needs some improvement. The problem is that X3D needs a built-in event-based mechanism to trigger complex, reusable animations. The widespread use of the event approach in web technology shows that it is an effective concept that adapts to different needs. The lack of such a mechanism in X3D makes creating animated and interactive scenes complex in some scenarios. References [1] P.X. Gao, A.R. Curtis,B. Wong,S. Keshav: It’s not easy being green. ACM SIGCOMM Computer Communication Review, 42(4), 211-22, 2012 [2] M.C. Hao, R.K. Sharma, D.A. Keim, U. Dayal, C. Patel, R. Vennelakanti: Application of Visual Analytics for Thermal State Management in Large Data Centres, Computer Graphics Forum, 29(6):1895-904, 2010 [3] S. Yan, X. Li: Analytical expression of indoor temperature distribution in generally ventilated room with arbitrary boundary conditions, Energy and Buildings, Vol. 208, 2020 [4] H.Edtmayer, D. Brandl, T. Mach, E. Schlager, H. Gursch, M. Lugmair, C. Hochenauer: Modelling virtual sensors for real-time indoor comfort control, Journal of Building Engineering, Vol. 67, 2023 [5] C. Lou, H.C. Zhou: Deduction of the two-dimensional distribution of temperature in a cross section of a boiler furnace from images of flame radiation, Combust. Flame, Vol. 143 (1–2), p.p. 97-105, 2005 JET 55 Gorazd Hren JET Volume 16 (2023) p.p. Issue 2, 2023 [6] G.Jiang, M.Kang, Z.Cai, H.Wang, Y.Liu, W.Wang: Online reconstruction of 3D temperature field fused with POD-based reduced order approach and sparse sensor data, International Journal of Thermal Sciences, Vol. 175, 2022 [7] B. Lange, N. Rodriguez, W. Puech, H. Rey, X. Vasques: A 3D particle visualization system for temperature management. Proceedings of SPIE - The International Society for Optical Engineering, 7868, 2011 [8] H.Edtmayer, D.Brandl, T.Mach, E.Schlager, H.Gursch, M.Lugmair, C.Hochenauer: Modelling virtual sensors for real-time indoor comfort control, Journal of Building Engineering. 67, 2023 [9] R. Jedermann, J. Palafox-Albarrán, P. Barreiro, L. Ruiz-García, J. Ignacio Robla, W. Lang, “Interpolation of spatial temperature profiles by sensor networks,” SENSORS, 2011 IEEE, Ireland, 778-781, 2011 [10] I. Cohen, D. Gordon: VS: a surface-based system for topological analysis, quantization and visualization of voxel data. Med Image Anal. 13(2), 245-56, 2009 [11] J. Kruger, P. Kipfer, P. Kondratieva,RD. Westermann: A particle system for interactive visualization of 3D flows, IEEE Trans Vis Comput Graph, 11(6), 744-56, 2005 [12] S.M. Seitz, C.R. Dyer: Photorealistic scene reconstruction by voxel coloring. Int J Comput Vision, 35(2):151-73, 1999 [13] N.A. Taranukha, Z.A. Izabekov: A method for voxel visualization of 3D objects. Program Comput Soft+; 33(6):336-42, 2007 [14] MatLab help. [15] R. Hajovsky, B. Filipova, M. Pies and S. Ozana: Using Matlab for thermal processes modeling and prediction at mining dumps, 12th International Conference on Control, Automation and Systems, p.p. 584-587, 2012 [16] Tim VOXview (https://www.mathworks.com/matlabcentral/fileexchange/78745-voxview), MATLAB Central File Exchange. Retrieved June 1, 2023 [17] F. Quintella, L.P. Soares, A.B. Raposo: DWeb3D: a toolkit for developing X3D applications in a simplified environment, Web3D ‘10: Proceedings of the 15th International Conference on Web 3D Technology, pp. 45–54, 2010 [18] X. Shen,H. Chen, T.M. Shih, X.Qingyu, Z. Hualin: Construction of three-dimensional temperature distribution using a network of ultrasonic transducers, Scientific Reports 9, 12726, 2019 [19] A. Schweigkofler, O. Braholli, S. Akro, D. Siegele, P. Penna, C. Marcher, L. Tagliabue, D. Matt: Digital Twin as energy management tool through IoT and BIM data integration, CLIMA 2022 conference, 2022 56 JET JET Volume 16 (2023) p.p. 56-68 Issue 2, 2023 Type of article: 1.01 http://www.fe.um.si/si/jet.htm SEAWATER AS A REAGENT IN THE FLUE GAS DESULPHURISATION PROCESS MORSKA VODA KOT REAGENT V PROCESU REAŽVEPLJEVANJA DIMNIH PLINOV Martin BriclR2 Keywords: Sea Water, Flue Gas Desulphurisation, Cleaning of Flue Gases, Sulphur Dioxide, Industrial Coastal Areas, Emissions Abstract Sulphur dioxide is a poisonous substance that is vastly formed in the process of fossil fuel burning inside the steam boiler of a thermal power plant, or any other industrial plant that uses fossil fuels as a prime source of energy. The technology involved in cleaning the sulphurous component from the raw, uncleaned flue gas flow in fossil fuel-based power and industrial plants is a form of technology that has been present for the last three decades and is constantly evolving in its characteristics and performance to deliver the highly-efficient flue gas cleaning procedure. The standard and technological mature technical solutions for flue gas desulphurisation of untreated flue gases are comprised of dry, semi-dry, and wet flue gas desulphurisation processes. The most frequently used solution, as well as that most applicable to most existing and new fossil fuelbased thermal power plants, is the wet flue gas desulphurisation process. The aforementioned wet cleaning process may be limestone-based (LFOS – Limestone Forced Oxidation System) or magnesium-based (MEL – Magnesium Enhanced Limestone). The paper focuses on the remaining wet flue gas desulphurisation procedure, one which is neither widely-present nor known within the industry – the SWFGD (Seawater Flue Gas Desulphurisation) process. This process has numerous advantages, including the presence of a costless reagent in vast amounts (seawater), the optimisation of the plant’s design, and those linked to operational costs. 2 Corresponding author: Dr Martin Bricl, mag.inž.str., Rudis d.o.o. Trbovlje, Trg revolucije 25b 1420 Trbovlje, Tel.: +386 3 56 12 409, E-mail address: martin.bricl@rudis.si JET 57 Martin Bricl JET Volume 16 (2023) p.p. Issue 2, 2023 We will present the flue gas cleaning process and its chemical aspects through a description of the process. We will then explore the main advantages and disadvantages of the corresponding process, as well as undertaking and presenting a comparative analysis between the main three wet flue gas desulphurisation processes (LFOS, MEL & SWFGD), taking into consideration all crucial points of each aforementioned wet flue gas cleaning process. As stated before, the main intention of wet desulphurisation processes is to remove the acid components from the untreated flue gas flow. In the process of doing this, the formation of by-products and effluent occurs, both of which have different impacts on the environment. Within the scope of this article, we will evaluate the environmental impact of the resulting by-products of each of the corresponding wet flue gas desulphurisation processes. Povzetek Žveplov dioksid je toksična snov, ki večinoma nastane pri izgorevanju fosilnih goriv v parnem kotlu termoelektrarne ali katerega koli drugega industrijskega obrata, ki uporablja fosilna goriva kot glavni vir energije. Tehnologija za čiščenje žveplove komponente iz neočiščenega toka dimnih plinov v termoelektrarnah in industrijskih obratih, ki delujejo na fosilna goriva, je prisotna že tri desetletja in se nenehno razvija v svojih lastnostih in zmogljivosti, da zagotovi učinkovit postopek čiščenja dimnih plinov. Standardne in tehnološko zrele tehnične rešitve za razžveplanje dimnih plinov so sestavljene iz suhih, polsuhih in mokrih postopkov razžveplanja dimnih plinov. Najpogostejši in najuporabnejši za večino obstoječih in novih termoelektrarn na fosilna goriva je postopek mokrega razžveplanja dimnih plinov. Zgoraj omenjeni postopek mokrega čiščenja je lahko na osnovi apnenca (LFOS – Limestone Forced Oxidation System) ali na osnovi magnezija (MEL – Magnesium Enhanced Limestone). V prispevku se bomo osredotočil na preostali postopek mokrega razžveplanja dimnih plinov, ki je v industriji premalo prisoten in poznan – postopek SWFGD (Sea Water Flue Gas Desulfurization). Ta postopek ima številne prednosti, kot so poceni reagent v velikih količinah (morska voda), optimizacija zasnove naprave in operativni stroški. Preko opisa postopka bomo predstavili potek čiščenja dimnih plinov in njegove kemijske vidike. Izpostavili bomo glavne prednosti in slabosti tega postopka ter izpeljali primerjalno analizo med glavnimi tremi postopki mokrega razžveplanja dimnih plinov (LFOS, MEL & SWFGD), ob upoštevanju vseh ključnih točk vsakega omenjenega mokrega postopka. Glavni namen razžveplanja dimnih plinov je odstraniti kislinske komponente iz neočiščenega toka dimnih plinov. Pri tem je prisotno nastajanje stranskih produktov in odplak, ki imajo različen vpliv na okolje, zato bomo v okviru prispevka tudi ovrednotili okoliške vplive nastalih stranskih produktov posameznega postopka mokrega razžveplanja dimnih plinov. 1 INTRODUCTION Due to the constant growth of the world’s population and, consequently, the growing demand for electricity, more and more energy plants are being built for the purpose of producing electricity. Despite the awareness that the combustion of fossil fuels releases greenhouse gases into our atmosphere, ultimately leading to the phenomenon known as the greenhouse effect, fossil fuel thermal power plants remain one of the key pillars of electricity production in the electricity supply sector. In addition to greenhouse gases, thermal power plants also emit other forms of emissions into our atmosphere, including sulphur dioxide (SO2). In high concentrations, sulphur dioxide emissions can be harmful to human health, affecting our living environment and plants, as well as other specimens within the environment itself. Over the last 30 years, the technology 58 JET Seawater as a reagent in the flue gas desulphurisation process involved in cleaning and removing acidic components from flue gases has advanced to a point where it can offer industrial and energy plants three ways to clean flue gases, namely through dry, semi-dry and wet flue gas cleaning processes. Each of these processes has its advantages and disadvantages, and the key to choosing the right flue gas cleaning technology is to carefully consider all factors that affect the technical, economic, and environmental parameters of the operation of the flue gas desulphurisation plant. This article covers and presents a basic record of wet flue gas cleaning technology using seawater as a reagent. [1] Such flue-gas desulphurisation plants are particularly suitable for locations where industrial or thermal power plants have access to large quantities of seawater. [2] 2 KEY MARKETS OVERVIEW Flue gas cleaning technology was introduced to the market in as early as 1970. The SWFGD process was first handed over to a Japanese client for commercial use in 1978 at an industrial plant for the production of chemicals. This was followed by a second SWFGD plant in 1988 at an oil refinery in Norway. [3] Despite good flue gas cleaning results and highly competitive investment and operating costs compared to other processes, this branch of flue gas cleaning technology has not gained a strong foothold in the market. [4] In 2000, the U.S. Environmental Protection Agency published an analysis of the use of FGD technologies. Figure 1: Areas where FGD technology is still in demand due to retrofit or green-field projects The results of the analysis showed that the share of SWFGDs represented only 0.6% of all desulphurisation systems in industrial and thermal power plants. Recently, we have seen an increase in demand and the installation of new SWFGD plants, mainly in coastal areas of Europe, Asia, and the Middle East (Turkey, Saudi Arabia, Iran, Iraq, Kuwait…). [5] Emphasis in the development of flue gas cleaning technology with the use of seawater mainly centres on reducing the space required for the installation of the plant, extending the life of the plant, and improving its design and operation to achieve better efficiency and reduce investment and operating costs. [6] JET 59 Martin Bricl 3 JET Volume 16 (2023) p.p. Issue 2, 2023 PROCESS DESCRIPTION The SWFGD plant consists of two main components. The first part is the system for the absorption of sulphur dioxide (SO2) and is intended to reduce the emissions of sulphur dioxide in the exhaust flue gases from the industrial or thermal power plant in question. The second part of the system is the process for the handling and treatment of used seawater before its discharge back into the sea. The principle of operation of SWFGD within a thermal power plant is schematically presented in Figure 2. Flue gases are channelled from the exit of the steam boiler to the electrostatic filter, the latter of which has the task of removing solid dust particles from the flue gases. Once the flue gases have passed through the electrostatic precipitator, they are further routed to the heat exchanger before entering the absorption vessel. In the heat exchanger, heat is transferred from hot, uncleaned, incoming flue gases to cold, cleaned outgoing flue gases from the absorption vessel. The cooled, uncleaned flue gases from the heat exchanger pass through into the absorption vessel. Here, the actual flue gas cleaning process also takes place. Figure 2: Principle of operation of the SWFGD within the thermal power plant In the absorption vessel, the flue gases from the steam boiler and the reagent - seawater come into contact. Most often, seawater is sprayed counter-currently through the flue gases. This allows optimal contact between the acidic components in the flue gases and the basic components in the seawater. The absorption zone in the absorption vessel consists of perforated steel plates that forcibly conduct flue gases so that contact between the aqueous and gaseous medium is maintained for as long as possible, consequently rendering the process more efficient. Cleaned flue gases are also conducted through a droplet separator that dries the flue gases. From the absorption vessel, the cleaned flue gases are directed into a heat exchanger, where they are heated to the appropriate temperature and then released through a chimney into the atmosphere. The reagent of the flue gas cleaning process - seawater - is directed to the top of the absorption vessel, where it is sprayed onto the flue gases. In doing so, a reaction takes place that separates sulphur dioxide [SO2] and hydrochloric acid [HCl] from the flue gases. Part or all of the reagent used in the purification process may be the cooling water originating from the condenser 60 JET Seawater as a reagent in the flue gas desulphurisation process in the Clausius-Rankin cycle. When the reagent reaches the bottom of the absorption vessel, it − + O→ HSO SO Hinto → spent HSO3−reagent + H + handling pool before being released back 2 + H 2and 3 2++H 2O the isSO retained then pumped (3.1) into the sea. When the flue gases come into contact with seawater inside the absorption vessel, − 2− +− − + begins +− the sulphur dioxide dissolve and + bisulfide is formed. Part of the resulting HSO →O SO + compound HHSO →O SO32− to +H 3 H 3 HSO 3H SO → H HSO +H (3.2) 2 + 2 then 3 2++ 2 to → 3This bisulfide can be SO converted sulphite. reaction occurs in the two chemical reactions (3.1) described below: 2− + − 2− + − + SO − + HSO+3− H →O SO + HSO HSO → 3 HSO 3 H 3 +H SO → 2 2 3 2++H 2O → HSO3 + H 2− + − + HSO3− → SO32− + HHSO 3 → SO3 + H (3.2) (3.1) (3.2) (3.2) Oxygen− is found both in flow of flue gases and+ in the seawater. Due to the presence of 2− the HSO3 + 12 O SO HSO +3−H+ +12 O2and →sulphite SO42− +from H the above two equations are oxidised to 2 → 4 bisulfide this oxygen, the resulting (3.3) sulphate through the following two chemical reactions: 2− − 1 2−2−2− − 1+ 2−2− + 1 SO SO 3 3++2 2OO 2 2→ 4 4 3 +3+ 44 + H HSO →SO SO HSO H+2 12OO2 2→ →SO SO (3.4) (3.3) 2− 2− 2− 2− − 1 2− − 1+1O → SO 2− + SO SOSO 3 ++2 1OO 2 → 4 3++ 4 HSO HSO +2 2 O2 2 → SO 3 2 → SO 4 3H 4 +H 2 (3.4) (3.4) (3.3) In the process of decomposition of sulphur dioxide in seawater to form bisulfide and sulphite − 2− 2− positively 2− 1 (and sulphate), charged hydrogen ions are also formed. These positively SO32consequently + 12 O2 → SO SO 4 3 + 2 O2 → SO4 (3.4) charged hydrogen ions oxidise the seawater in the reaction vessel and lower its pH. Thus, acidified seawater must be properly neutralised to maintain the appropriate conditions for the existence and preservation of ecosystems and organisms in seawater. Neutralisation of previously used seawater inside the absorption vessel is achieved with compounds that form an integral part of seawater. The chemical reactions describing the neutralisation of acidified seawater are as follows: − HCO3− + H + → CO HCO OH + → CO2 + H 2O 2 +H 3 2+ (3.5) (3.6) (3.6) (3.5) To achieve the appropriate alkalinity, the neutralisation process takes place in the first zone of + − 2− + − 2− − H +→ − H seawater +→ HCO the collection where fresh from the cooling system of the thermal power plant CO CO 3 ++ 3 3 ++ 3 HCO H basin, →HCO CO HCO (3.6) 3 2 +H 3 2OH → CO2 + H 2O (3.5) is added to the oxidised seawater. The medium in the first part of the collection basin must be thoroughly mixed before entering the next zone. A further part of the collection basin is the zone − 2− + − where into the well-mixed seawater medium. Blowing air and thoroughly mixing CO32− the + Hair+ is→blown HCO CO 3 3 + H → HCO3 (3.6) the effluent from the SWFGD process constitutes an important step from an ecological and technical point of view if we want to ensure the effective operation of the treatment process. [7] By optimally blowing air into the pool, we achieve: − 2− − − 2− − CO HH+ +→→HCO CO +O HH+ +→→HCO 3 3++ 3 H HCO CO HCO CO23 + H 2O 23 + 3 2+ • Optimal and satisfactory oxidation and compounds; • Adequate efficiency of neutralisation of the medium before its discharge into the sea; • The restoration of adequate oxygen concentrations in seawater. JET 61 Martin Bricl 4 JET Volume 16 (2023) p.p. Issue 2, 2023 SEAWATER AS THE PROCESS REAGENT Seawater used as a reagent in the flue gas cleaning process is subject to certain requirements before being discharged back into the sea. The following are particularly important: • When the process seawater is taken from the sea itself, the pH value of the seawater returned from the process back to the sea must be as close as possible to the original pH value; • The temperature of the seawater being returned to the sea must not significantly exceed the average seawater temperature before its initial removal from the sea; • The COD (Chemical Oxygen Demand) value must be as low as possible; • The value of dissolved oxygen must be suitably high; • The level of sulphates in the seawater being returned to the sea must be as low as possible; • The concentration of particulate matter in the seawater being returned to the sea must be correspondingly low. Upon the discharge of the absorption vessel, once the flue gas cleaning has been completed, the pH value of the reagent – seawater - is about 3 to 4. Since the critical value for the acidity of seawater for marine and marine organisms is about 6.5, acidified seawater must be deacidified before its discharge into the sea. The basicity of seawater is increased through the aforementioned process of blowing air into the seawater before releasing it back into the sea. At the outlet from the pool, where the air is forced into the used seawater, the pH value is controlled, which must be between 6 and 7. When seawater returns from this process, it is about 2 °C warmer prior to its discharge into the sea compared to when it was initially taken from the sea. The amount of dissolved oxygen in seawater capture before entering the process is between 50% and 100%, while at the outlet, following its treatment in the mixing basin and after the forced air blowing process, the proportion of dissolved oxygen in seawater is closer to 70% to 90%, or more than 6 mg/L. Near thermal power plants using the seawater desulphurisation process, long-term seawater wastewater quality monitoring has not demonstrated any adverse effects on the ecosystem, the marine environment, or on the animals. The so-called COD (Chemical Oxygen Demand) method is also used to determine the quality or pollution of seawater at the outlet of the mixing and ventilation basin. This is the standard method for indirectly determining the degree of contamination of an aqueous sample. This method is based on the chemical decomposition of organic and inorganic compounds and elements in an aqueous sample. The result of this method gives us the equivalent of the amount of water - dissolved oxygen in the aqueous sample (given in ppm) or in mg/L, which is used for its decomposition by water pollutants during a two-hour decomposition process in boiling potassium dichromate. In the case of a smaller industrial plant with an SWFGD plant, the permitted number of COD in the discharged seawater is approximately 100 - 150 mg/L. If the SWFGD plant is upgraded to a larger thermal power plant (700 MWe and more), where the flow of consumed seawater as a reagent through the flue gas cleaning system can exceed 100,000 m3/h, which is equivalent to the flow of a smaller river, the COD limit is 5 mg/L. 62 JET Seawater as a reagent in the flue gas desulphurisation process Table 1: Typical chemical elements & compounds within seawater Dissolved oxygen in water allows the oxidation of sulphite to sulphate. The rate of this reaction is influenced by the pH value of the seawater. The most efficient flue gas cleaning procedure inside the absorption vessel should be at a pH value of seawater between 4.1 and 4.5 and between 5 and 5.6. Most of this oxidation takes place inside the absorption vessel. However, for complete oxidation, it is necessary to forcibly blow oxygen into the sea-wastewater treatment pool. With the SWFGD plant, it is possible to achieve high efficiency of sulphur extraction from the flue gas stream in thermal power plants (> 99%), but it is also necessary to ensure that the COD value (Chemical Oxygen Demand) is between 2.5 and 5.0 mg/L. The oxidation process continues even after the discharge of processed seawater back into the sea, and so it is important to ensure an adequate level of oxygen upon the discharge of the medium into the sea. [8] The amount of oxygen in the processed water discharged into the sea during the SWFGD process is regulated by Council Directive 79/923/EEC. In Table 2, the discharged seawater quality parameters of a typical SWFGD process are presented. Parameter Unit Inlet seawater Discharge / 8.2 6-7 COD mg/L O2 0 2.5 - 5.0 DO % 50 - 100 70 - 90 mg/l 2700 2785 Temperature ⁰C T T + 1,5 Salinity % 33.5 33.5 mg/L SS SS + 1 pH Sulphate Suspended solids The following section describes and discusses the degree of impact stemming from the main factors that may provoke environmental concerns, as well as their potential consequences on relevant organisms. The rise of temperature in seawater, where effluent water is returned to the sea, occur in small area where temperature is elevated in duration less than 1 minute. Therefore, we can assume that any rise in the water’s temperature will have an insignificant JET 63 Martin Bricl JET Volume 16 (2023) p.p. Issue 2, 2023 impact on the environment, the quality of the seawater, and the species residing in it across the broader observed area. However, the rise in local sea temperature at the moment of discharge of the effluent water back into the sea represents the primary factor of mortality for the pelagic organisms (see Table 3). The pH – pH of effluent water returned to the sea can have a pH between 6 to 7. Pelagic organisms are exposed to this pH for a period of 15 minutes, depending on the average current velocity. This time exposure to pH in the range of 6 – 7 should not cause any significant effect on pelagic organisms. Benthic organisms will also not experience any adverse effects, since the effluent water from the SWFGD process is returned to the sea and dispersed in the upper layers of the seawater. Chemical Oxygen Demand and dissolved oxygen – the expected values of both parameters are presented in Table 2. Organisms can be affected by a deficit of dissolved oxygen within 250 m from the discharge of effluent water from the SWFGD process. The direct exposure to dissolved oxygen deficit is evaluated for 15 minutes. Salinity – in the case of seawater evaporation in cooling towers and scrubbers, the 12% increase in the concentration of Na+ ions can be expected in the effluent water entering the collection pit and aeration channel. Many of the pelagic organisms can adapt to variation in the salinity of the seawater in the range of ± 18%. Therefore, the possible change in salinity concentration should not have a direct impact on pelagic organisms. Sulphate – ambient seawater contains an amount of SO4 in the range of 3200 – 3650 mg/L, depending on the amount of salinity concentration in the local seawater. As stated in the description of salinity change, a variation in sulphate concentration is also acceptable for the pelagic and benthic organisms within the range of ± 18%. Suspended solids – the desulphurisation process using seawater as a reagent produces an additional amount of suspended solids in the effluent water returned to the sea within a range of less than 1mg/L. Table 3: Survival possibility of species in nearby seawater where SWFGD scrubber effluent is discharged back to the sea* Species Milkfish larvae Milkfish juveniles Copepods Clams Eel elvers Black seabreams Tiger prawns * Effluent concentration □ □ □ □ □ □ Seawater temperature □ □ □ □ □ □ Legend: □ – not affected; - affected None of the aforementioned environmental impacts originating from the SWFGD process are severely harmful to the environment, nor the species that reside within it. We can thus conclude that seawater, when used as the medium in the scrubbing process of raw flue gas flow in the SWFGD process, is a good medium for scrubbing sulphur components in raw flue gases before they enter our atmosphere. [9] Seawater is especially appropriate to be used as a reagent since it has an inherited alkalinity concentration of approximately 100 – 110 mg/L of CaCO3. 64 JET Seawater as a reagent in the flue gas desulphurisation process At the same time, it is important to point out that seawater already has an existing amount of sulphur in it, approximately 0.9 kg of sulphur per 1 tonne of seawater. [10] 5 CONSTRUCTION MATERIALS The choice of construction materials is crucial to ensure a suitably long service life and to safeguard the initial investment amount of the flue gas treatment plant. High corrosion in flue gas desulphurisation plants can be caused by both the flue gases originating from the steam boiler and the oxidised seawater that is retained in the absorption vessel. The use of Hastelloy® and Duplex stainless steel is very expensive and can greatly increase the initial investment in an SWFGD treatment plant. To optimise investment costs, other construction materials can also be used, such as FRP (Fibre Reinforced Polymer), chlorine-butyl rubber, silicone rubber, etc. It is also possible to use ordinary structural steels, which are wall-coated with stainless steel. It is crucial to foresee the proper material choice for key process equipment. Only by doing this can the longlasting and reliable operation of flue gas desulphurisation plants be achieved. The main absorber vessel and pumps with corresponding pipelines constitute key pieces of equipment. 6 COST COMPARISON OF ABATEMENT TECHNOLOGIES The following technologies for cleaning flue gases in a big thermal power plant in terms of removing sulphur dioxide from raw flue gases are different in their operating principles as well as in initial investment and, later, their operational costs. Table 4: Cost comparison for the flue gas desulphurisation process using different available abatement technologies DESULPHURISATION OF FLUE GASES WET FGD WET FGD UPGRADE SEMI DRY FGD DRY FGD SWFGD Load Factor [%] Emission Rate [mg/Nm3] BREF limit [mg/Nm3] Capital cost [€/MW] 30 70 30 < 100 < 100 100 150-175 150-175 130-320 300.000,00 300.000,00 22.500,00 Cost per unit power generation [€/ MWh] 7,00 € 5,99 € 2,13 € 70 100 130-320 22.500,00 3,90 € 30 100-500 130-320 90.000,00 3,62 € 70 100-500 130-320 90.000,00 5,11 € 30 70 30 70 50-200 50-200 < 200 < 200 130-320 130-320 130-320 130-320 110.000,00 110.000,00 70.000,00 70.000,00 3,37 € 3,88 € 0,83 € 1,98 € JET 65 Martin Bricl JET Volume 16 (2023) p.p. Issue 2, 2023 Figure 3: Cost of abatement technologies The new BREF directive is going to limit the emission rates for big thermal power plants as shown in table 2. Depending on the steam boiler load factor, the new BREF limits will need to be met. In the case of the new wet FGD system, the new regulative proposes emission limits between 150 – 175 mg/Nm3. In the case of upgrading the existing wet FGD system, the emission limits propose the emissions in the range between 130 – 320 mg/Nm3. The same emissions level range is predicted for the semi-dry and dry process of flue gas desulphurisation processes as well as for the SWFGD. [11] 7 CONCLUSION Since the SWFGD process necessitates access to fresh seawater, it is locally limited to coastal areas. This same flue gas cleaning process can also be used in the use or incineration of coal with a high sulphur content, but consideration must be given to the fact that high concentrations of sulphur dioxide in the flue gases would require an additional amount of seawater in the process to achieve adequate efficiency. This additional required amount of seawater cannot be provided only from the cooling system of the thermal power plant, and so it is necessary to provide additional captures of seawater directly from the sea, something which increases the investment and operating costs of flue gas treatment plants. The main advantages of the SWFGD process are as follows: • SWFGD does not typically require any chemical reagents for its operation; • During the process of operation of the treatment plant, no by-products are created that could be dangerous for the environment and living beings; • The design of the plant is simple. This enables the plant’s reliable operation and long service life, assuming that the appropriate construction materials were selected during its construction; 66 JET Seawater as a reagent in the flue gas desulphurisation process • The process allows up to 99% removal of sulphur dioxide from flue gases. For adequate economic profitability of the SWFGD process, we envisage the use of fuel - coal, with a sulphur content of up to 1.5%. If coal with a higher sulphur content is burned in a thermal power or industrial plant, then the additional use of sorbents is required within the SWFGD flue gas cleaning system. These sorbents can be limestone, sodium hydroxide, or magnesium oxide. The use of such sorbents is sometimes necessary to achieve adequate removal of sulphur dioxide from the flue gases in the absorption vessel and to neutralise seawater waste before it is discharged from the process and is returned back to the sea. References [1] Oikawa, K.: Seawater Flue Gas Desulphurization: Its Technical Implications and Performance Results, Environmental Progress, Vol. 22, No. 1, April 2003 [2] Zhang, Y.: Effects of Additives on Seawater Flue Gas Desulphurization, International Conference on Environment Science and Engineering, Vol. 8, Singapore 2011 [3] Abrams, J. Z.: Use of Seawater in Flue Gas Desulphurization, JAPCA, Vol. 38, No. 7, 1988 [4] Quian, Z.: Non-calcium desulphurization technologies, Clean Coal Centre, Vol. 170, June 2010 [5] Ushiku, T.: MHPS’s Environmental Technology in Asian market, Mitsubishi Hitachi Power Systems Ltd., Vol. 220, No. 8401, February 2014 [6] Sargent & Lundy: Wet Flue Gas Desulphurization Technology Evaluation, National Lime Association, January 2003 [7] Bakke, T.: Evaluation of the environmental impact from flue gas seawater scrubbing at the Qatalum aluminium smelter in Qatar, Report of the Norwegian Institute for Water Research, September 2010 [8] Lee Twu, J.: Using Seawater to Remove SO2 in a FGD System, Waste Water - Treatment and Reutilization, 2011 [9] Maikandaan, T. P.: Development of Flue Gas Desulphurization Reactor for Reducing SOX Emission from Flue Gases Emitted in Thermal Power Plants, International Journal of Innovative and Exploring Engineering (IJITEE), Vol. 9, No. 2S4, December 2019 [10] Guo, H.: A Numerical Investigation on the Optimization of Uneven Flow in a Marine De-Sox Scrubber, Processes, Vol. 8, No. 862, July 2020 [11] Watanabe, Y.: Development and Installation of Marine-use Hybrid Sox Scrubber System that Complies with IMO Sox Emission Regulations, Mitsubishi Heavy Industries Technical Review, Vol. 53, No. 2, June 2016 JET 67 Martin Bricl Nomenclature (Symbols) (Symbol meaning) LFSO Limestone Forced Oxidation System MEL Magnesium Enhanced Limestone SWFGD Sea Water Flue Gas Desulphurisation SO2 Sulphur Dioxide U.S. United States FGD Flue Gas Desulphurisation HCl Hydrogen Chloride COD Chemical Oxygen Demand 68 JET JET Volume 16 (2023) p.p. Issue 2, 2023 Author instructions Author instructions http://www.fe.um.si/sl/jet.html http://www.fe.um.si/si/jet.htm MAIN TITLE OF THE PAPER SLOVENIAN TITLE Author1, Author 2, Corresponding author  Keywords: (Up to 10 keywords) Abstract Abstract should be up to 500 words long, with no pictures, photos, equations, tables, only text. Povzetek (Abstract in Slovenian language) Submission of Manuscripts: All manuscripts must be submitted in English by e-mail to the editorial office at jet@um.si to ensure fast processing. Instructions for authors are also available online at http://www.fe.um.si/en/jet/author-instructions.html. Preparation of manuscripts: Manuscripts must be typed in English in prescribed journal form (MS Word editor). A MS Word template is available at the Journal Home page. A title page consists of the main title in the English and Slovenian language; the author(s) name(s) as well as the address, affiliation, E-mail address, telephone and fax numbers of author(s). Corresponding author must be indicated. Main title: should be centred and written with capital letters (ARIAL bold 18 pt), in first paragraph in English language, in second paragraph in Slovenian language. Key words: A list of 3 up to 6 key words is essential for indexing purposes. (CALIBRI 10pt) Corresponding author: Title, Name and Surname, Organisation, Department, Address, Tel.: +XXX x xxx xxx, E-mail address: x.x@xxx.xx  1 Organisation, Department, Address 2 Organisation, Department, Address 78 JET JET 69 1 2 Authors names and surnames Paper title JET Vol. 16 (2023) Issue 1 ---------Abstract: Abstract should be up to 500 words long, with no pictures, photos, equations, tables, text only. Povzetek: - Abstract in Slovenian language. Main text should be structured logically in chapters, sections and sub-sections. Type of letters is Calibri, 10pt, full justified. Units and abbreviations: Required are SI units. Abbreviations must be given in text when first mentioned. Proofreading: The proof will be send by e-mail to the corresponding author in MS Word’s Track changes function. Corresponding author is required to make their proof corrections with accepting or rejecting the tracked changes in document and answer all open comments of proof reader. The corresponding author is responsible to introduce corrections of data in the paper. The Editors are not responsible for damage or loss of submitted text. Contributors are advised to keep copies of their texts, illustrations and all other materials. The statements, opinions and data contained in this publication are solely those of the individual authors and not of the publisher and the Editors. Neither the publisher nor the Editors can accept any legal responsibility for errors that could appear during the process. Copyright: Submissions of a publication article implies transfer of the copyright from the author(s) to the publisher upon acceptance of the paper. Accepted papers become the permanent property of “Journal of Energy Technology”. All articles published in this journal are protected by copyright, which covers the exclusive rights to reproduce and distribute the article as well as all translation rights. No material can be published without written permission of the publisher. Chapter examples: 1 MAIN CHAPTER (Arial bold, 12pt, after paragraph 6pt space) 1.1 Section (Arial bold, 11pt, after paragraph 6pt space) 1.1.1 Sub-section (Arial bold, 10pt, after paragraph 6pt space) Example of Equation (lined 2 cm from left margin, equation number in normal brackets (section.equation number), lined right margin, paragraph space 6pt before in after line): Equation 70 JET (1.1) JET 79 3 Paper title Authorsnames names and surnames Authors and surnames ---------- JET Volume 16 (2023) JET Volume 16 (2023) p.p. Issue2, 1, 2023 2023 Issue Tables should have a legend that includes the title of the table at the top of the table. Each table should be cited in the text. Table legend example: Table 1: Name of the table (centred, on top of the table) Figures and images should be labelled sequentially numbered (Arabic numbers) and cited in the text – Fig.1 or Figure 1. The legend should be below the image, picture, photo or drawing. Figure legend example: Figure 1: Name of the figure (centred, on bottom of figure, photo, or drawing) References [1] N. Surname: Title, Journal Title, Vol., Iss., p.p., Year of Publication [2] N. Surname: Title, Publisher, Year of Publication [3] N. Surname: Title [online], Publisher or Journal Title, Vol., Iss., p.p., Year of Publication. Available: website (date accessed) Examples: [1] J. Usenik: Mathematical model of the power supply system control, Journal of Energy Technology, Vol. 2, Iss. 3, p.p. 29 – 46, 2009 [2] J. J. DiStefano, A.R. Stubberud, I. J. Williams: Theory and Problems of Feedback and Control Systems, McGraw-Hill Book Company, 1987 [3] T. Žagar, L. Kegel: Preparation of National programme for SF and RW management taking into account the possible future evolution of ERDO [online], Journal of Energy Technology, Vol. 9, Iss. 1, p.p. 39 – 50, 2016. Available: http://www.fe.um.si/images/jet /Volume 9_Issue1/03-JET_marec_2016-PREPARATION_OF_NATIONAL.pdf (7. 10. 2016) Example of reference-1 citation: In text [1], text continue. Nomenclature (Symbols) t 80 JET (Symbol meaning) time JET 71 72 JET 9 771855 574008 JE T I Jo ur na l o f E ne r g y Te chno lo gy I Vo l. 16, I s s ue 2, Se p te mb e r 2023 I Unive rsit y o f M a rib o r, Fa cu l t y o f E n e rgy Te c h n o l o gy ISSN 1855-5748