**Supply Chain Management for Bioenergy and Bioresources**

Special Issue Editors **Dionysis Bochtis Charisios Achillas**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Special Issue Editors* Dionysis Bochtis Institute for Bio-Economy and Agri-Technology (iBO) Greece

Charisios Achillas Technical Educational Institute of Central Macedonia Greece

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Energies* (ISSN 1996-1073) (available at: https://www.mdpi.com/journal/energies/special issues/ Bioenergy Biorecourses).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03943-462-6 (Pbk) ISBN 978-3-03943-463-3 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Special Issue Editors**

**Dionysis Bochtis** works on the area of systems engineering focused on bio-production and related provision systems with enhanced ICT and automation technologies up to fully robotized production systems. Throughout his carrier he has held various positions, including Professor in Agri-Robotics at the Lincoln Institute for Agri-Food Technologies, University of Lincoln, UK, and Senior Scientist in Operations Management at the Department of Engineering of Aarhus University, Denmark. Currently, he is the Director of the Institute for Bio-Economy and Agri-Technology (IBO), Centre for Research & Technology Hellas (CERTH). He runs numerous research projects on ICT and robotics in agricultural production. He is author of more than 300 articles (100 in peer-reviewed journals) and has been invited for more than 30 keynote speeches around the globe. He is involved in the boards of a number of international scientific organizations, holding positions including the following:


**Charisios Achillas** is Assistant Professor at the International Hellenic University, Department of Supply Chain Management. He graduated in 1999 from the Department of Engineering (Aristotle University of Thessaloniki), with a degree in Mechanical Engineering. He continued his studies with an MSc in Engineering Project Management in 2001 (UMIST, Manchester, UK). In 2009, he received his PhD Doctorate in the field of reverse logistics. Dr. Achillas is a senior researcher at the Institute of Bio-Economy and Agri-Technology, Center for Research & Technology Hellas (CERTH) and the Laboratory of Heat Transfer and Environmental Engineering, Department of Mechanical Engineering, Aristotle University of Thessaloniki. Since 2003, Dr. Achillas has been involved in science, research and development, from technical development to coordination, including project and financial management. His work has flourished in participating in more than 40 research projects (the majority being EU-funded), dealing mostly with environmental engineering, reverse logistics and sustainable development. Dr. Achillas is the author of more than 170 scientific publications. His research interests span the fields of environmental management, sustainable development and circular economy, with a primary focus on reverse logistics and agri-chains.

## **Preface to "Supply Chain Management for Bioenergy and Bioresources"**

In the modern world, the competitiveness of bioenergy- and/or bioresources-related activities heavily depends on the effectiveness of supply chain management. A large number of multidisciplinary topics are involved in the bioresources and bioenergy production fields. Although the technical issues that are related with the topic are well-discussed and do not represent major barriers, supply chain management issues, such as design of the network, collection, storage or transportation of bioresources, are still considered as fundamental questions that need to be answered to enable the optimal exploitation of bioenergy and bioresources. Moreover, modeling of material and energy flows; identification of the dynamic character of the supply chains; available reverse logistics (waste management) alternatives; economic, social and environmental sustainability of bioresource supply chains; novelty in the applied business models; and decision support frameworks towards efficient supply chain management for bioenergy and bioresources present critical operational sustainability issues and business-making potential.

This Special Issue, entitled "Supply Chain Management for Bioenergy and Bioresources", includes one extensive review on yellow and woody biomass supply-chain management, together with six original papers which span around a number of innovative, multifaceted, technical developments that are related to all different echelons of supply chain management for bioenergy and bioresources. More specifically, in their work ("Green, Yellow, and Woody Biomass Supply-Chain Management: A Review"), Rodias et al. present a comprehensive review on research studies targeting biomass supply-chain management advancements related to three types of biomass sources, namely green biomass sources (such as perennial grasses), yellow biomass sources (such as crop residues) and woody biomass sources (such as willow). From their review, it becomes evident that the presented up-to-date trends on biomass supply-chain management and the potential for future advanced application approaches play a crucial role in business and sustainability efficiency of a biomass supply chain.

The Special Issue is also enriched by the following six studies: (1) Lampridi et al. ("Energy Footprint of Mechanized Agricultural Operations") present a modeling methodology for the precise calculation of the energy cost of performing an agricultural operation, with their model incorporating operational management into the calculation, while simultaneously considering the commercially available machinery (implements and tractors); (2) Vlachokostas et al. ("Decision Support System to Implement Units of Alternative Biowaste Treatment for Producing Bioenergy and Boosting Local Bioeconomy") propose a generic methodological scheme based on multicriteria analysis for the development of small-, medium- or large-scale units of alternative biowaste treatment, with an emphasis on the production of bioenergy and other bioproducts, taking into account environmental, economic and social criteria to support robust decision-making; (3) Papageorgiou et al. ("Forecasting of Day-Ahead Natural Gas Consumption Demand in Greece Using Adaptive Neuro-Fuzzy Inference System") examine the application of neuro-fuzzy models, so as to develop a real, accurate natural gas (NG) prediction model for Greece, thus providing a fast and efficient tool for utterly accurate predictions of future short-term natural gas demand; (4) Vamvanoulas et al. ("Dry Above Ground Biomass for a Soybean Crop Using an Empirical Model in Greece") propose a new empirical equation for the estimation of daily dry above ground biomass for a hybrid of soybean, which presents a useful tool for estimations without using destructive sampling; (5) Papageorgiou et al. ("Decision-Making Process for Photovoltaic Solar Energy Sector Development using Fuzzy Cognitive Map Technique") focus on the investigation of certain factors and their influence on the development of Brazilian photovoltaic solar energy with the help of fuzzy cognitive maps, an established methodology for scenario analysis and management in diverse domains, inheriting the advancements of fuzzy logic and neural networks; (6) Wu et al. ("A Cloud-Based In-Field Fleet Coordination System for Multiple Operations") analyze the structure and composition of an auto-steering-based collaborative operating system for fleet management that is developed to realize an in-field flow-shop working mode often adopted by the scaled agricultural machinery cooperatives.

In conclusion, this Special Issue seeks to contribute to the bioenergy and bioresources agenda through enhanced scientific and multidisciplinary knowledge that may boost the performance efficiency of supply chain management and support the decision-making process of stakeholders. We are confident that the papers included in the present Special Issue provide interesting case studies that may motivate researchers and encourage additional efforts in supply chain management for bioenergy and bioresources.

> **Dionysis Bochtis, Charisios Achillas** *Special Issue Editors*

## *Review* **Green, Yellow, and Woody Biomass Supply-Chain Management: A Review**

**Efthymios Rodias 1,2, Remigio Berruto 1,\*, Dionysis Bochtis 2, Alessandro Sopegno <sup>1</sup> and Patrizia Busato <sup>1</sup>**


Received: 8 July 2019; Accepted: 1 August 2019; Published: 6 August 2019

**Abstract:** Various sources of biomass contribute significantly in energy production globally given a series of constraints in its primary production. Green biomass sources (such as perennial grasses), yellow biomass sources (such as crop residues), and woody biomass sources (such as willow) represent the three pillars in biomass production by crops. In this paper, we conducted a comprehensive review on research studies targeted to advancements at biomass supply-chain management in connection to these three types of biomass sources. A framework that classifies the works in problem-based and methodology-based approaches was followed. Results show the use of modern technological means and tools in current management-related problems. From the review, it is evident that the presented up-to-date trends on biomass supply-chain management and the potential for future advanced approach applications play a crucial role on business and sustainability efficiency of biomass supply chain.

**Keywords:** supply-chain design; strategic planning; operational planning; energy crop production; crop residue

#### **1. Introduction**

For the generation of any product, a sequence of processes connected to design, decision making, and execution, and a series of financial, information, and material flows are performed throughout different stages of the production. The different stages of the production constitute an integrated system called supply chain. The basic aim for the successful design of a supply chain is to meet the requirements of the final customer regarding a specific product. The supply-chain concept for agri-food and biomass relates to not only manufacturing and retailing sectors but also to agricultural sector.

In recent decades, energy crops constitute a highly-potential share among crops taking into account the need for greener energy production. Energy crops are crops that are cultivated for biomass, biogas, or other biofuels (e.g., biodiesel, bioethanol) production. They are mostly green crops that come from wild nature, such as perennial grasses with high potential for bioenergy production. Green-type biomass includes crops such as *Miscanthus*, *Panicum virgatum* (also known as switchgrass), *Arundo donax*, etc. At the same level, yellow biomass refers to crop residues that come from any crop and represent another category of biomass production related to feedstock. Examples of yellow biomass are corn stover, wheat straw, etc. It should be mentioned that the main sources of biogas production are the energy crops and the use of agricultural residues [1]. On top of this, there is a variety of woody crops that contribute significantly to biomass energy production globally. Woody biomass is any biomass that is connected to wood sources. Examples of woody biomass sources are willow, poplar

short-rotation coppice, etc. All of them have various and complex constraints regarding the entire management policy and practices that should be followed for the optimal biomass production.

The main challenges regarding supply chain issues on each biomass type category, as they are presented above, are different. Challenges related to green-type biomass include any grass-type crop operational issues, including particular issues in harvesting and handling (such as optimal scheduling), and less on other processes (such as soil cultivation or fertilization) due to their easy adaptability to various environments. On the other hand, yellow-type biomass requirements include optimal collection, handling and transportation processes. This incorporates on-time scheduling of collection and transportation and optimal task execution in cases where multiple fields are covered. Challenges on woody biomass sources are different from the other two types given its operational processes particularities. The woody energy crops require different crop establishment, cultivation, harvesting, and transportation processes. Of course, some operational issues would be similar with the other two types (such as scheduling of operations), but there are technical issues that are solved in different ways. A short example regarding collection and transportation would be about harvesting of green or yellow biomass in bale form compared with woody crops that whole trees are collected.

At this time, supply-chain management (SCM) in agricultural production, handling, and transportation processes is vital and there are always various issues that should be faced through better SCM. There is a large amount of research works regarding SCM of green, yellow, and woody biomasses. The work here, targets on an up-to-date literature review on recent publications on green, yellow, and woody biomass SCM.

The main challenges for the creation of this review were, firstly, to underline variation practically in different approaches for specific existing problems in SCM, secondly, to provoke the development of supply-chain management on the specific target group by proposing possible solutions on various upcoming matters and, finally, to provide a brief review of various followed practices/methodologies and their effects on the SCM.

There are previous reviews on the supply-chain management in agricultural processes regarding green, yellow, and woody biomass types. A biomass supply chain evaluation and optimization are suggested by a literature review regarding forest feedstock [2]. A systematic review was presented in order to present the key factors throughout the biomass supply chain of green and woody crops that affect the application of bioenergy buffers in complex bioenergy production systems [3]. A wider review about biofuel SCM is conducted under the objective of uncertainties and sustainability issues [4]. On the opposite side, a more practical review was presented regarding many types of mathematical models in bioenergy crops production, including both energy crops production processes and transportation but also biorefinery/biomass conversion modeling processes [5]. Even though the scientific contribution of these reviews is highly important, to the knowledge of the authors, there is no recent review regarding chain management aspects for green, yellow, and woody supply.

The objective of this paper is to highlight and focus on green, yellow, and woody biomass supply-chain management research works (52 studies), and, as a second step, to create a classification in order to propose opportunities for further research by focusing on research gaps and identified issues needed to be tackled.

#### **2. Supply-Chain Management Definition**

In order to conclude to a successful definition of SCM, i is important to make a short comparison between traditional management (TM) practices and SCM practices. Under the financial concept, by TM a reduction in a company's costs may be achieved, while by SCM a whole-chain cost efficacy will be obtained. Regarding data exchange and monitoring information, the first case is limited on the business's own needs, while in SCM it can be extended for whole-chain planning and/or monitoring processes. Another point is the coordination between different levels of a channel, where in TM there is only a single contact for the interchange among the channel pairs, while by SCM multiple contacts and coordination between various businesses and levels of channels can be accomplished. Finally, there is a number of risks and rewards that cannot be shared under the TM philosophy, but by SCM all the risks and rewards are shared in a long-term period. The above-mentioned comparison is only a small sample of the differences between TM and SCM, suggesting the need for assimilating more and more SCM practices.

There are various ways to describe and define SCM. For the purpose of this paper a definition of the term regarding crop production processes would be: Supply-chain management is the integrated planning of in-field and/or logistics operations, application of these operations, coordination between the different levels of the channel(s) and, finally, control of all processes and necessary activities in order to produce and transport, in the most efficient way, the products that finally will satisfy the requirements of a given market [6].

Given this definition, we could set that all the in-field and logistics processes are included in the term supply chain, not only as a physical operation but also as a decision-making activity, both associated by material flows and exchange data and, as a consequence, the correlated financial/energy flows. In this light, the supply chain includes not only the producer and its suppliers, but also includes the processing units, logistics operators, warehouses, etc.

#### **3. Review Methodology**

The methodology followed in this review includes a series of theoretical considerations taken into account in the pre-processing stage. Throughout our analysis, the included steps (as also presented in Figure 1) are the following:


#### *3.1. Eligibility Criteria*

Studies related to green, yellow, woody, and multiple-type biomass SCM are included in this study. At the same time, as woody biomass is only referred to woody energy crops (such as willow), publications related to forest species or forest waste biomass are excluded from the scope of this review. The SCM-related tasks are included in the wide boundary in which this review's scope is subjected, as presented in Figure 2. For the scope of this study, the included publications should be related to SCM, and should be referred mainly to the crop production processes and/or the transportation from farm to the field and/or from field to the plant. Any further biomass conversion operations or processing was considered to be out of the scope of this review. However, publications that refer to individual operations and not to inter-connected parts of the supply chain were not taken into account.

A number of literature-related eligibility criteria are set. These criteria include: (1) The work should be published in English language, (2) the included studies should be research articles published in peer-review journals, and (3) they should have been published within the last five years (current one included; i.e., from 2015 up to present). Publications in journals that are not research articles, such as reviews or short communications are excluded from this review paper.

**Figure 1.** The summary of the review process.

**Figure 2.** Processes included in this review's boundary.

#### *3.2. Information Sources*

There are a total of 52 studies included, which withdrawn from the electronic databases: Web Of Science, ScienceDirect, Wiley Online Library, and SpringerLink. The primary searches were implemented in the mid-end of January 2019. It is possible that any newly emerging publication after this date is not included. Nonetheless, the keywords, their combinations, and the searching strings are saved in order to make the whole process replicable for future use.

In the present review, no grey literature is included (i.e., any literature produced in electronic or print format that has not been controlled by any commercial publisher, for example, technical or research reports, doctoral dissertations, etc.).

#### **4. Biomass-Type Classification Frameworks and Analysis**

The existing literature regarding green, yellow, or woody biomass sources was classified in the following categories: (1) Literature focusing on supply-chain strategic planning, and (2) literature regarding operational planning for a series of operations (or the whole system) throughout the supply chain. In the term "supply-chain strategic planning", system productivity or life cycle analysis (LCA) (or energy consumption/balance assessment) of the biomass crop under study was considered. In a similar way, any financial/economics evaluation of a certain crop production/operational processes would be included under supply-chain strategic planning. However, these categories are not rigid and there is high possibility a particular work could be incorporated in both of them. A brief summary of the literature allocated to these categories according to the specific biomass type reported is presented below.

#### *4.1. Green Biomass*

#### 4.1.1. Strategic Planning

Studies related to green biomass feedstock that include assessment of the system productivity or LCA analysis (including sub-categories such as energy consumption/balance, etc.) or financial evaluation are incorporated in this category. Apart from this, studies that include innovative design strategies throughout supply chain are also under consideration here.

An up-to-date geographical information system (GIS)-based approach by using a supply-chain simulation model is presented [7]. The objective of the presented approach was the determination of the optimal locations of beet crops that maximizes the profit of bioethanol production plants. Another study on allocation of energy crops to dispersed field locations is the one presented by [8], where a crop-to-field allocation tool targeting the maximum energy gain of the system was proposed. In that study, energy consumption is considered for all in-field operations and transportation in the routes between farm-field and field-plant. The optimization problem was modeled both as a linear and as binary programming. Another simulation approach focuses on spatial geographical allocation of *Panicum virgatum* crop fields, storage facilities, and biorefineries were presented [9]. The work was based on a two-phase modeling process under the scope of economical sustainability of the supply chain components. This two-phase simulation included the implementation of the agricultural land management alternative with numerical assessment criteria (ALMANAC) model for crop productivity simulation combined with the AnyLogic simulation model that is capable, among others, of finding the optimal locations for biomass storage facilities.

In [10], alternative supply chain configuration scenarios in *Miscanthus* harvesting and transport were compared, financially optimized, and assessed from a sustainability perspective for different annual demand in biomass, yield, and time of harvesting, among others. For the same energy crop, a computational tool was presented by [11], based on an in-depth analysis and energy requirements estimation of individual *Miscanthus* fields, including all the in-field operations and transportation. An economic approach based on the implementation of GIS was implemented to evaluate green biomass (*Miscanthus* and *Panicum virgatum*) SCM under different structures of the supply chain [12]. The scale of examination was downsized from

county-level to field-level in order to more easily find the suitable individual areas for biomass production and, in this way, allowing microeconomic evaluation and modeling.

*Miscanthus* was also under the scope of study in terms of supply-chain strategic planning considering the production operations [13]. A simulation model called Simulateur mulTIdiscplinaire pour les Cultures Standard (STICS) was presented (as an improvement of an existing one) for the accurate prediction of the produced biomass in the long-term (up to 20 years from planting time) in various case studies. On top of this, the authors evaluated the nitrogen content in different cropping environments and in-field management practices. An approach based on soil maps, climatic data, and observed yields of *Miscanthus* fields was presented by [14] and in combination with GIS data concluded in the perspective of a decision support tool development for optimal supply-chain strategic planning. An optimization tool was also evaluated in the *Miscanthus* production supply chain and presented promising results as a decision support tool for farmers—in crop cultivation strategy development—and for policy makers—in monitoring and improvement of supporting practices [14].

A combined financial and energy requirements analysis was conducted by [15] for harvesting and transportation operations in *Arundo donax* production systems. Targeting the optimal strategic planning of supply chain, different operational alternatives were proposed.

Sugarcane is evaluated in [16] study, which focuses on sugarcane SCM and also on green harvesting residues management under different operational practices. The authors developed a simulation model to present the biomass flow through the whole supply chain (in-field and transport operations). Except for operational constraints, they also took into account weather and geographical constraints to identify possible bottlenecks in the supply chain and biomass availability, such as non-synchronized harvesting, handling, and transporting operations.

Triticale is a hybrid of wheat and rye not widely considered as a crop for bioenergy production. A study that focused on the assessment of triticale as a potential biorefinery feedstock is the one presented by [17]. They introduced improved harvesting methods that have a positive effect on costs reduction and availability secure of high-quality biomass in the long-run. A number of autochthonous perennial grasses, as potential crops for biomass production, were evaluated in terms of supply chain efficiency [18]. Authors conducted a four-year experiment to examine parameters such as energy efficiency, crop yield, and water-use efficiency.

Energy balance estimation comes to the center of interest for complex biomass production systems, as presented by [19]. More analytically, they presented an assessment process for the energy balance of multiple-crops and multiple-fields systems by implementing a web-based tool. Three different crops were evaluated as a case study, namely corn silo, wheat, and rapeseed. In [20], multiple-production systems were also evaluated by developing and implementing a comparative computational tool potentially applied to any set of given crops, given energy-related or production-related input, and according to any specific production practices. A study that refers to different bioenergy cropping systems is presented by [21]. They evaluated the productivity of selected systems by using experimental data that were further analyzed by a simulation tool to identify potential limitations of resource-use efficiency on biomass dry matter yield.

Carbon footprint assessment throughout the supply chain is included also in supply-chain strategic planning. An integrated assessment of commercial crops (such as *Maize*) with bioenergy crop production (such as switchgrass) can be an innovative way for biomass supply-chain strategic planning, as presented by [22]. Authors took into consideration the effect of SCM practices on cost, yield production, and system sustainability. A multiple-crop CO2 annual emissions estimation regarding cultivation processes was conducted by [23]. In addition to this, they made a long-term (for a 30-year period) prediction of the CO2 emissions for the selected crops in order to identify the most sustainable crop(s), from the environmental point of view, for bioethanol production. Environmental impact assessment of multiple perennial green biomass crops in marginal land was presented by [24], based on crop comparison under specific conditions and further development and improvement of the supply chain.

#### 4.1.2. Operational Planning

A multiple-optimization strategy was followed considering the impacts of operational management on minimization of costs, maximum, yields, and sustainability of *Panicum virgatum* supply chain [22]. For this scope, authors modeled the integration of *Panicum virgatum* on a real corn production field and its effects on profitability, productivity, and environmental improvements of the system by using mainly the landscape environmental assessment framework (LEAF) tool (Version 2.0, United States).

*Arundo donax* is also the target energy crop in study [25], where authors provided an agronomical assessment testing on the crop structure and regrowth potential as they are affected by the harvesting time. They further compared different harvesting alternatives and evaluated the biomass quality based on the harvesting parameters, such as harvesting time and harvesting frequency.

Sweet sorghum and sugar beet are two important energy crops for bioethanol production. In this light, many studies have focused on the SCM of these crops. The biomass yield of sweet sorghum and sugar beet was estimated in an experimental study in southern Italy [26]. In the same work, the energy performance under different in-field operational management scenarios was evaluated for these two crops. The effect of three levels of shade (low, medium, and high) on *Maize* production (such as growth and yield) that is cultivated for biogas production was evaluated in [27]. These levels of artificial shading were presented to affect the biogas-related parameters (such as leaf area index and energy availability for plant growth) and the final biogas and methane yield. In a *Miscanthus* production system, biomass production, costs, and supply-chain constraints are considered in a whole supply chain financial optimization strategy [10].

There is also reference considering specific aspects of supply-chain operational issues in green biomass crops. The optimized design of the in-field operations following optimal route planning (B-patterns) under the objective of minimizing energy cost was assessed in two green cropping biomass production systems (*Miscanthus* and *Panicum virgatum*) [28].

#### *4.2. Yellow Biomass*

#### 4.2.1. Strategic Planning

The establishment of biogas plants in the optimal location is directly connected to crop residue-related parameters (such as quantity, accessibility, weather conditions, etc.). Regarding the strategic planning of yellow biomass supply chain, the models multiple linear regression and artificial neural networks (ANNs) were implemented for the estimation of available crop residues (specifically, corn stover and wheat straw) in multiple sites [29]. In the same work, potential suitable locations for bioenergy plants are identified by this approach, and subsequently, the optimal plant location is suggested together with the biomass delivered cost. For the same crop residues, authors in [30] calculated the potential sustainable yellow biomass quantities while maintaining soil productivity and health. They focused on large-scale bioenergy applications and proposed the establishment of sustainable removal rates of residues and supply chain cost for various regions.

A simulation-based model that considers multiple-locations assessment and selects the optimal location for bioethanol plants based on parameters such as wheat biomass density, supply chain network, etc., was presented by [31]. The focus was on the financial analysis and environmental impact assessment.

Regarding crop residues from corn production (stover), different biomass handling and transportation scenarios were evaluated under an LCA analysis in [32]. The selected scenarios included combinations of biomass handling (in bale form or pelletized form) and either storage in an intermediate depot or by directly transportation from field to the biorefinery.

Cotton stalks represent another type of crop residue biomass. In [33] an integrated GIS and an ANN high spatial resolution model was developed to assess available cotton stalks harvesting and transportation. In addition, by GIS analysis the suitable biorefinery locations were selected under the criterion of the minimum total transport distance and the delivered cost. Finally, the estimation of spatial and temporal variations of potential cotton stalks in the United States was presented by [34], including also an assessment of different SCM practices for cotton stalks under an engineering economics approach.

A GIS-based model taking into account various features connected to a geographic region in order to determine the location of biomass-processing facilities considering also uncertainty factors [35]. These factors are connected to the availability of crop residues, other environmental, financial, weather, or even social constraints. In this study, residues coming from various crops are used, such as corn, sorghum, sugarcane, wheat, barley, agave, rice, and pecan nut. The development of a spatial decision tool is also the main objective in [36], with the aim of advancing this research field. The developed tool is useful at the strategy development level where the decision making is conducted. The optimal alternative solutions or combination of them were determined by identifying the best regions for the establishment of biomass plants, according to the crop residue availability, by including also factors such as geographical features, accessibility of the crop residues points, and the complexity and structure of the logistics network. Multiple types of crop residue biomass were also under consideration in [37]. Authors presented a location-route collaborative optimization model that they used to design the multiple residues collection network and corresponding routes. The yellow biomass considered was stalk, leaves, etc., from crops such as cereals and beans, tobacco stalks and leaves, stalk from beet, sugarcane, among others. The sites were depicted by nodes and thus the node capacity of the network and the logistics cost were the basis for the main modeling structure. On top of that, they proposed and built a mathematical algorithm to solve the location-routing problem under specific constraints.

#### 4.2.2. Operational Planning

Fruit orchards, vineyards, and olive trees are of high biomass interest given the produced residues in the field. On this, [38] presented an integrated framework for the development of a decision support tool for logistics operations of fruit trees, vineyards, and olive tree pruning and branches. In this work, the necessity of a "smart logistics system" is presented. Additionally, the components of this system under certain technical and other requirements and constraints were defined together with the determination of the information to be managed by the system. The olive tree-pruning biomass supply chain was accurately estimated and optimized on supply-chain development and modeling frameworks [39].

Innovative techniques that include quantification of pruning biomass residues coming from individual olive trees and based on LiDAR (light detection and ranging) data processing methods were used by [39]. Their work is of high value for supply-chain operational planning given that they focused on the accurate estimation of the biomass volume and distribution and other aspects such as supply-chain modeling and development.

#### *4.3. Woody Biomass*

#### 4.3.1. Strategic Planning

SCM is quite significant also for woody biomass considering strategic planning. Authors in [40] analyzed the viability of eucalyptus supply chain and further evaluated the optimal location of a power plant to be fed by wood. Their analysis was based on GIS databases to calculate available biomass under certain restrictions. Even though this work is one among many that target the optimal bioenergy plant localization, their methodology includes calculation of optimized costs and CO2 emissions throughout the supply chain considering local and regional data. Authors in [41] developed a discrete event simulation model for techno-economic assessment of harvesting, handling, and transport operations for the supply of short-rotation coppice willow (SRCW) to energy conversion plants. Their main objective was to suggest cost-effective supply chain solutions in order to deliver all-year-round SRCW to the biomass plant.

Poplar short-rotation coppice (SRC) belongs also to woody species under consideration for biomass production. In the work presented in [42], the most important environmental impacts of wood chips production from poplar SRC on marginal land were provided in detail, and a comparison was made between common SCM practices and alternative SCM practices in terms of modified harvesting and handling systems, irrigation, and fertilization.

Eucalyptus is included in the woody biomass SCM category and is considered one of the most commercial woody crops cultivated for bioenergy production. In [43], the economic viability of various operational practices regarding whole eucalyptus tree trunk production was analyzed under varying production and economical parameters. To analyze risk management a Monte Carlo analysis was run. Stochastic models also included different planting densities and arrangements in order to provide a framework for assessment of financial risk applied in decision analysis for bioenergy crop investments.

Willow was under consideration in study [44], where part of the supply chain was modeled and evaluated (harvesting, collecting, and transportation). Authors run the integrated biomass supply analysis and logistics (IBSAL) simulation model (Oak Ridge National Laboratory, Oak Ridge, TN, United States) in a significant number of SR willow fields and evaluated the effect of the major input factors (i.e., parcel size, field shape, crop yield, field-to-plant distance, and type of collection machinery) in the performance of the system.

An innovative mobile application for application in mostly woody biomass supply-chain strategic planning was suggested by [45]. The objective was to create a useful tool for the farmers in order to be able to evaluate the energy potential of tree pruning in their own orchards, share the information regarding the specific region that the orchard belongs, and match this data with the availability for heating energy in agro-industrial or other buildings.

#### 4.3.2. Operational Planning

In [40], optimization of supply chain was pointed out, targeting the optimal location of wood-fired plants under specific constraints. SRC willow, as it has been presented above in [41], is under the scope of supply-chain operational planning (regarding harvesting and handling operations) in the supply chain from field to plant. Additionally, the impact of the most significant input factors was evaluated regarding the total willow supply chain performance [44].

Poplar tree production was also evaluated on part of the supply chain related to cultivation and harvesting operations [46]. The main objective was to point out the environmental impact of poplar trees cultivated in various soil nutrient levels and suggest the optimal strategy. Similarly, in poplar SRC, optimization of cultivation and transportation practices through a financial and environmental analysis was presented by [47].

Regarding woodchip transportation on a short supply chain (up to 70 km distance), an analysis of the energy input and the carbon dioxide emissions was conducted by using two types of transportation means (i.e., agricultural tractors or industrial trucks) and carrying out a thorough analysis under parameters related to road design (such as, traffic lights, intersections, etc.), traffic conditions, and road surface conditions (such as rain, fog, ice) [48].

#### *4.4. Multiple Biomass*

#### 4.4.1. Strategic Planning

Several studies include combination of green and yellow biomass supply-chain strategical and operational planning. A work that combines green and yellow biomass SCM based on agro-ecosystem models was presented by [49]. The main objective of the work was the supply chain assessment and design (including crop production and transportation) under the demand of a bioethanol plant with feedstock from three types of crops grown for different purposes, namely: crop residues; annual crops; and perennial crops (i.e., wheat straw, triticale, sorghum, and *Miscanthus*). The developed model focused on productivity measures and the environmental impact of the supply chain. Among others, the scope of the model was the minimization of the balance between food and feed crops. Under the scope of farmers' profit maximization and of environmental sustainability, a landscape approach for the

assessment of multiple various-sized subfields was carried out in [50]. They evaluated opportunities coming from yellow biomass collection and transportation together with green biomass potential development on low-perspective fields that normally make no profit.

Supply-chain strategic planning is underlined by energy balance estimation of in-field and transportation operations research studies regarding green and woody biomass types, such as the ones presented by [19] and [51].

A comprehensive model that can be applied for strategic planning of green, yellow, and woody biomass supply chain was presented in [52], where authors developed a mixed-integer linear programming model that targeted the minimization of the establishment cost of the integrated supply chain under specific constraints (e.g., biomass availability, capacity of technology, etc.). Authors divided the supply chain in two parts: The biomass supply chain that referred to the part from field to biorefinery, and the liquid fuel supply chain that referred to flows of ethanol and blended fuels from biorefinery to gas station. The focus of the developed model was on supply-chain-related solutions and biomass-related solutions that may contribute to a more financially viable supply chain.

Finally, a CyberGIS platform (CyberGIS Gateway, Illinois, United States) was developed by [53] as a decision support tool for biomass supply-chain management and analysis. More specifically, authors described, in their model, how to integrate optimization within the CyberGIS platform in order to solve complex spatial decision-making problems throughout the supply chain.

#### 4.4.2. Operational Planning

Another study combined green and yellow biomass assessment [54]. Authors evaluated both production and transportation operations of corn (stover) and *Panicum virgatum* under two different baling systems, and, in a second step, they compared the corresponding supply chains financially, in terms of energy cost, and in terms of GHG emissions.

Operational management practices were assessed within three different designed supply chains of green biomass in terms of operational times and costs [55]. In particular, three different approaches where simulated and the effect of yield uncertainty and machinery productivity were evaluated based on a sensitivity analysis approach. Authors in [56] proposed a simulation approach for the optimal biomass bales location. Various bale aggregation methods were simulated and evaluated in terms of the required in-field transport operations for different methods. Similarly, on [57] authors simulated the in-field layout of biomass bales, providing the desired bale stacks layout under various different parameters (such as field area and shape, yield, bale mass, etc.) for increasing the efficiency of biomass handling and transportation operations.

Scheduling of handling and transport operations is considered crucial in any biomass supply-chain strategic planning. In this light, a tool was developed for scheduling of field machinery used in these operations in geographically-dispersed fields under specific constraints [58]. In that work, authors presented an algorithm that creates individual working schedules for multiple machines that carry out more than on consecutive operations at multiple fields, taking into account the given readiness in each specific field on the specific timing. They used a basic Tabu Search method and a modified one that produces optimized work planning.

#### **5. Methodological Approach-Based Classification Framework**

The research studies, which were already classified above in two levels, can be further classified according to the methodological approach followed in each cited work. In this light, the four categories of included studies (i.e., green, yellow, woody, and multiple-type biomass) were kept the same as in the previous classification framework. The studies related to green-type biomass are summarized and presented in Table 1.



*Energies* **2019**, *12*, 3020

**Table 1.** *Cont.*

Another significant share of the potential biomass sources for bioenergy production regards yellow-type biomass. The studies that include agricultural residues SCM—only those that come from crop production and not any agricultural residue—are presented in brief in Table 2. These studies include various methods such as simulation and experimental methods, but also real-time artificial intelligence methods that are connected to machine learning technologies; technologies that have lately emerged in the agricultural sector [59].

Similarly, with the previous two literature categories, Table 3 regards studies correlated to woody biomass SCM approaches that come from woody crops but not common forest species. Moreover, in Table 4 the studies that refer to SCM of more than one different type of biomass feedstock are listed.

#### **6. Discussion and Conclusions**

In this review, 52 research studies were included in total, published in 24 different journals (Figure 3). "Biomass and Bioenergy" represent the most referenced journal in the current review followed by "BioEnergy Research" and "Journal of Cleaner Production".

**Figure 3.** Number of papers extracted by each journal.


studiesclassifiedbycropresidueandapproach.


Woody-typebiomasscitedstudiesclassifiedbycropand


Generally, the use of final biomass product may vary, either for heating or electricity production closer to its primary natural form (i.e., pure biomass), either by converting biomass to biofuels (such as bioethanol, biodiesel, or biogas) through biorefining processes. Here, three main categories (namely biomass, biogas, bioethanol) as types of final product were pointed out in the included literature and presented in Figure 4, classified further to each different type of biomass. According to the literature, in the case of multiple types of biomass sources, the final product mostly leads to direct biomass use and less to bioethanol production. In parallel, literature shows that green and yellow biomass types may be equally used for biogas production. Regarding the bioethanol-related works, all types of biomass are referenced except woody type. This is also representative of the existing situation in bioethanol and biogas production industries.

**Figure 4.** Biomass type in relation to the final product.

The reviewed literature was assessed under the methodological approach defined above and the results were presented in Tables 1–4 for each biomass-type. An aggregated representation of this categorization is presented in Figure 5. In general, various types of simulation approaches were used in almost 43% of the total studies of this review, while 26% used experimental methodological approaches and 31% used other types of approaches, such as online tools, web applications, etc.

As previously mentioned, the parts of the supply chain that were included in this review correspond to the stages from in-field production processes up to transportation to the biomass processing facilities, including any intermediate steps of the chain. However, the steps of the supply chain included in each conducted study varied. On this, eleven different combinations of supply chain steps were found in total in all biomass types. The number of works allocated to each one of these categories is presented in Figure 6. There is a significant amount of studies—19% of total—that correspond to in-field production and biomass transport stages, while 17% of the studies are targeted to in-field processes.

**Figure 6.** Number of works allocated to various parts of supply chain: (**1**) Harvest and transport; (**2**) production; (**3**) production and transport; (**4**) collection and quantification; (**5**) collection and transport; (**6**) handling and transport; (**7**) harvesting, collection, transport, and handling; (**8**) transport; (**9**) harvesting, collection, and transportation; (**10**) harvesting and handling; (**11**) harvesting, handling, and transport.

Among the articles, three supply chain categories are related to energy crops biomass production and six categories are connected to crop residues biomass supply-chain management. On top of this, there are seven categories related to woody energy crops biomass supply chain and six categories are presented to connect with studies that take into account more than one type of biomass (such as combined green- and yellow-type).

Overall, 52% of the referenced publications are related to strategic planning of supply chain, 29% of the works are related to operational planning, while 19% of the works present combined methods (targeting strategic and operational planning of supply chain). Green biomass-related literature includes approaches related to strategic planning (12 works), operational planning (five publications), while there are four works that refer to both strategic and operational planning tasks of supply chain. Regarding yellow biomass literature, five and one approaches are, respectively, targeted to strategic and operational planning, while three works include both strategical and operational planning Similarly, for woody biomass literature, three works are related to strategic planning, four are related to operational planning, and three are related to approaches that combine strategic and operational planning of biomass supply chain. Finally, in literature that includes multiple biomass types, strategic planning of supply chain is targeted seven times, five works refer to operational planning, and three works are related to a combination of strategic and operational planning approaches.

Regarding the methodology-based framework, as it was presented in Figure 5, simulation techniques are the mostly used methodological tool in green-type biomass and multiple-type biomass systems. In yellow and woody type biomass, experimental approaches and other methods (such as analytical models, etc.) are the main methods.

Combining simulation methodological approaches to experiment-based methods in any kind of biomass type may result in the real evolution of these systems under the scope of optimal SCM. Moreover, in the future, methodological tools and other types of approaches are expected to contribute to this evolution by developing innovative real-time tools that would be useful, especially in solving problems in complex biomass systems. There is already a trend in this direction but this is quite dynamic and can be widely expanded more in many levels, both for the needs of the mentioned crops and systems but also in other biomass systems. By this review, an integration of different approaches is suggested for further evolution in currently-used methodology under the scope of production levels increase, sustainability enhancement of the systems at hand, and for ensuring bio-based product quality.

**Author Contributions:** Writing—Original Draft Preparation, E.R., D.B. and R.B.; Methodology, E.R., D.B. and P.B.; Investigation, E.R., D.B. and A.S.; Conceptualization D.B. and R.B.; Writing—Review and Editing, E.R. and D.B.; Supervision, R.B. and D.B.

**Funding:** This research received no external funding.

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

#### **References**


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

## *Article* **Energy Footprint of Mechanized Agricultural Operations**

#### **Maria Lampridi 1,2,\*, Dimitrios Kateris 1, Claus Grøn Sørensen <sup>3</sup> and Dionysis Bochtis <sup>1</sup>**


Received: 19 December 2019; Accepted: 7 February 2020; Published: 10 February 2020

**Abstract:** The calculation of the energy cost of a cultivation is a determining factor in the overall assessment of agricultural sustainability. Most studies mainly examine the entire life cycle of the operation, considering reference values and reference databases for the determination of the machinery contribution to the overall energy balance. This study presents a modelling methodology for the precise calculation of the energy cost of performing an agricultural operation. The model incorporates operational management into the calculation, while simultaneously considering the commercially available machinery (implements and tractors). As a case study, the operation of tillage was used considering both primary and secondary tillage (moldboard plow and field cultivator, respectively). The results show the importance of including specific operation parameters and the available machinery as part of determining the accurate total energy consumption, even though the field size and available time do not have a significant effect.

**Keywords:** agricultural operations; energy use; assessment tool; workability; machinery

#### **1. Introduction**

The concept of agricultural sustainability has been a highly debated subject in recent years, not only because of its importance with respect to evaluating the agri-food system, but mostly because of its multivariate nature [1]. In fact, there is a large number of parameters that are considered with respect to the primary production of agricultural products [2]. These parameters include not only the relevant material flows (e.g., water, fuel, fertilizers, seeds, and crops), but also the means of production, which include labor and the use of agricultural machinery along with the inclusion of operation management [3]. The integration of the individual evaluation of each of those parameters with respect to their contribution to sustainability concerns more and more researchers in recent years. However, agricultural sustainability is approached in a variety of ways in the literature, examining, individually or in combination, economic, environmental, and social indicators [4]. These indicators attempt to assess the impacts of agriculture in a quantitative or a qualitative manner. Several indicators as well as indicator indexes and frameworks have been proposed in the literature with respect to each dimension of sustainability [5]. The most frequently used economic indicators include profitability, income, efficiency, productivity and yield, while the societal indicators include the education, health risk, employment, operational difficulties and access to health [6–9]. With respect to the environmental indicators, the most popular are energy use, water footprint, soil quality, biodiversity, and greenhouse gas (GHG) emissions [10–13].

In particular, energy cost (or energy consumption or energy input) is one of the most examined agricultural sustainability indicators of field operations [14]. Energy consumption is an important agricultural environmental indicator and serves as a basis for the calculation of several composite indicators related to energy such as energy efficiency and energy balance [14,15]. The research focus during recent years was on the assessment and calculation of the energy footprint for the entire crop production, providing approaches for improvement in terms of environmental impact measures, through the minimization of agricultural inputs such as water use, fertilizers, and plant protection products [14–17]. In all of the cases, the goal is to assess or minimize the environmental burden of agricultural operations, while maintaining the quality and quantity of production [13]. For these assessments, the calculation of the energy consumption was mostly based on reference values and databases, while the models used do not always consider the commercially available equipment and mostly focus on material consumption [18]. Additionally, the effect of machinery selection or machinery dimensioning (considering, for example, the time availability, that is, the workability of the system at hand) is not considered during the assessment, while there are only a few studies examining the embodied energy of agricultural machinery [19–21] and the energy cost of machinery when performing a specific task [22]. It is also worth noting that, with respect to the selection of the required machinery, mostly financial methods are employed that consider the economic cost of operations [23,24]. However, with the advancement of technology it is important to examine the energy contribution of the machinery, especially when the introduction of new equipment is assessed [25–27].

Considering the above and attempting to fill the identified research gap in the field of agricultural machinery energy assessment, the present study attempts to incorporate "machinery dimensioning notion" in the calculation of the energy costs. The ultimate goal is to assess the influence of the workability that determines the available time window for performing a specific field operation on the resulting energy consumption for this operation. Furthermore, the machinery used is assessed, examining the contribution of the standardized-commercial machinery in the overall energy requirements. For that reason, a computational model was designed and implemented in order to examine the machinery energy cost of performing an agricultural operation. The model selects the required machinery to perform the operation under given workability constraints [23], choosing from a collection of commercially available implements and tractors, and considering a series of operational parameters including the soil type, the operation depth, and the field area to be processed.

In the following sections, the methods and the detailed design of the model are presented followed by results from two case studies. The examined operations concern the soil preparation and, more specifically, light soil cultivation with a field cultivator and primary tillage with a moldboard plow.

#### **2. Materials and Methods**

#### *2.1. Energy Cost of Agricultural Operations*

The execution of an agricultural operation requires both direct and indirect inputs. To that end, the employed approach in this study borrows from the well-established methodology of financial cost calculation [28], and allocates the total energy required to perform an operation to direct and indirect energy (Figure 1). The indirect energy concerns the energy assessment of inputs that are not related to the operation performed but concern the acquisition and the ownership of the equipment. The direct energy involves the energy assessment of the direct inputs related to agricultural operations.

On the basis of the above, the indirect energy is differentiated between the machinery and repair and maintenance embodied energy (MJ · kg<sup>−</sup>1), as well as the housing embodied energy (MJ · <sup>m</sup>−<sup>2</sup> · <sup>y</sup>−1). The machinery embodied energy refers to the energy required to manufacture the agricultural equipment [20,29]. The embodied energy of repair and maintenance can be either expressed as a percentage of the manufacturing energy requirements [30,31] or can be estimated in absolute numbers, as attempted by Mantoam et al. [32]. The housing embodied energy represents the energy contribution of all the building and infrastructure necessary for the storage of the equipment [30]. It should be noted that these denoted parameters are usually calculated for the entire life cycle of the machinery or building [33]. In this study, the energy contribution from manufacturing, repair and maintenance, and housing is reduced to the duration of the operation performed considering as related to the total life of the equipment.

**Figure 1.** Energy costs of agricultural operations.

The operation energy is divided into the fuel and lubricant energy (MJ · l −1 ), and the embodied energy of labor (MJ · <sup>h</sup>−<sup>1</sup> ). The energy of fuels and lubricants is related to the fuel and lubricant consumption as calculated based on the ASABE (American Society of Agricultural and Biological Engineer) standards [34]. The embodied energy of labor expresses the energy consumed by the workers performing the operation [35]. Several approaches have been proposed in the literature regarding the estimation of the embodied energy of labor differing in the system boundaries and operations considered [16,36–39].

#### *2.2. Energy Estimation Model*

In order to calculate the energy cost of performing an agricultural operation, the model presented in Figure 2 was used. The calculation consists of two distinct stages, namely the machinery selection and the calculation of the energy costs. For the machinery selection stage, in order to estimate the effect of standardized machinery on the overall energy consumption, two models were considered, namely the continuous and the discrete model, while the energy cost calculation is the same for both models. The continuous model calculates the energy consumption of a theoretical machine system without correction to standard commercial machinery, while the discrete model considers the commercially available implements and tractors for the determination of the energy costs. More specifically, the necessary machine system (single or multiple if it is required) to perform the operation is automatically calculated by the model, considering both operational parameters and the machinery available. The calculation process of the discrete model is presented in detail below, however, the basic equations are standard for both models. Nevertheless, it should be highlighted that, in the continuous model, only a single machine system is considered.

During the machinery selection stage, first the available time tav (h) to perform an agricultural operation is calculated (Equation (1)):

tav = twp · twt · Wcoef (1)

where *twp* (*d*) is the working period in days, twt - <sup>h</sup> · <sup>d</sup>−<sup>1</sup> is the working time, and Wcoef (−) is the workability coefficient [40]. Given a specific available time, the theoretical area capacity Cth - <sup>m</sup><sup>2</sup> · <sup>h</sup>−<sup>1</sup> is calculated (Equation (2)):

$$\mathbf{C\_{th}} = \frac{\mathbf{A}}{\mathbf{t\_{av}}} \tag{2}$$

where A m2 is the size of the field under examination. The theoretical capacity determines, as a sequence, the theoretical minimum width wth (m) (Equation (3)) that must be applied in order to complete the operation on time.

$$\mathbf{w}\_{\text{th}} = \frac{\mathbf{C}\_{\text{th}}}{\mathbf{E}\_{\text{f}} \cdot \mathbf{s}} \tag{3}$$

where Ef (%) stands for the efficiency of the operation performed, while s - km · <sup>h</sup>−<sup>1</sup> is the operating speed.

**Figure 2.** Calculation model.

Then, according to the discrete model, the theoretical width is adjusted based on the commercially available implements. When the theoretical width is smaller than the maximum width available, then the width selected is the one that is immediately bigger than the theoretical width. When the theoretical width is larger than the maximum width available, it means that more than one implements are required to perform the operation on time. In that case, the maximum width is selected (as many times as required), while the remaining width is adjusted to the next bigger available width. In should be noted that, in this step, the actual number of machines required to perform the operation is determined. Considering the above, the actual width wact (m) of the operation performed is formulated and, as a consequence, the actual area capacity is calculated by Equation (4) as follows:

$$\mathbf{C\_{act}} = \mathbf{w\_{act}} \cdot \mathbf{s} \cdot \mathbf{E\_{f}} \tag{4}$$

Finally, the actual duration of the operation is calculated (Equation (5)):

$$\mathbf{t}\_{\text{act}} = \frac{\mathbf{A}}{\mathbf{C}\_{\text{act}}} \tag{5}$$

In the case of the continuous model where a correction to standard machines is not performed, the actual width and actual capacity are the same as the theoretical width and capacity and actual time coincides with the time available to perform the operation. Once the required implements have been determined, the calculation of the required tractors follows. The implement draft Ddraft (N) is calculated according to the ASABE standards [41] by Equation (6):

$$\mathbf{D}\_{\rm draft} = \mathbf{F}\_{\rm i} \cdot \left[ \mathbf{A} + \mathbf{B} \cdot \mathbf{s} + \mathbf{C} \cdot \mathbf{s}^2 \right] \cdot \mathbf{w}\_{\rm act} \cdot \mathbf{O}\_{\rm depth} \tag{6}$$

where Fi (−) is the soil texture adjustment parameter; A, B, and C are machine specific parameters; and Odepth (cm) is the operation depth. Considering the above, the draft power (kW) is calculated as follows [23]:

$$P\_{\text{draft}} = \frac{D\_{\text{draft}} \cdot s}{1000} \tag{7}$$

Finally, the total power required Prequired (kW) is calculated from Equation (8) as follows:

$$P\_{\text{required}} = \frac{P\_{\text{draft}}}{\text{Tr}\_{\text{coef}} \cdot 0.83} \tag{8}$$

where Trcoef (−) stands for the tractive coefficient [42], while the coefficient 0.83 is used for the Power Take Off (PTO) to gross-flywheel conversion [34]. The calculated Prequired (kW) corresponds to the minimum power required to perform the operation considering the necessary implements as they were defined in the previous steps. However, commercial tractors are available with a standard power capacity for each model, as defined by the manufacturer. To that end, the required power is adjusted to the next bigger available. In the case of the continuous model, the required power coincides with the actual power because there is not a correction to standard tractors. The machinery selection stage provides the list of tractors-implements that are necessary for the completion of the procedure in the given time available for the discrete model and the single tractor and implement system for the continuous model. For each tractor-implement set, the energy costs of performing the operation are calculated.

Regarding the operation parameters calculations, the hourly fuel consumption Fcon <sup>l</sup> · <sup>h</sup>−<sup>1</sup> is calculated from Equation (9) (for diesel) as follows [34]:

$$\mathbf{F}\_{\rm con} = \left( 2.64 \cdot \mathbf{P}\_{\rm ratio} + 3.91 - 0.203 \cdot \sqrt{738 \cdot \mathbf{P}\_{\rm ratio} + 173} \right) \cdot \mathbf{P}\_{\rm required} \tag{9}$$

where Pratio expresses the ratio of equivalent PTO power necessary for an operation to that maximum available from the PTO [34]. The relevant lubricant consumption Lcon- <sup>l</sup> · <sup>h</sup>−<sup>1</sup> is calculated in Equation (10) (for diesel) [34].

$$\text{L}\_{\text{con}} = 0.00059 \cdot \text{P}\_{\text{required}} + 0.02169 \tag{10}$$

Lastly, it is important to note that, for the calculation of the operation energy costs, the transportation duration is also considered along with the operation duration. The resulting per hectare total energy EnTot (MJ) is calculated as the sum of the direct EnDirect and indirect energy EnIndirect, as presented in Equation (11):

$$\text{En}\_{\text{Tot}} = \text{En}\_{\text{Direct}} + \text{En}\_{\text{Indirect}} \tag{11}$$

The direct energy EnDirect (MJ) is calculated by Equation (12) as follows:

$$\text{En}\_{\text{Direct}} = \text{En}\_{\text{Fuel}} + \text{En}\_{\text{Lurbr}} + \text{En}\_{\text{Labor}} \tag{12}$$

where EnFuel (MJ) is the total fuel energy; EnLubr is the total lubricant energy and EnLabor is the total labor energy expressed in Equations (13)–(15):

$$\mathbf{E}\mathbf{n}\_{\text{Fuel}} = \mathbf{F}\_{\text{com}} \cdot \mathbf{t}\_{\text{act}} \cdot \mathbf{E}\mathbf{n} \mathbf{b}\_{\text{Fuel}} \tag{13}$$

$$\text{En}\_{\text{Lurbr}} = \text{L}\_{\text{con}} \cdot \text{t}\_{\text{act}} \cdot \text{Emb}\_{\text{Lurbr}} \tag{14}$$

$$\text{En}\_{\text{Labor}} = \text{t}\_{\text{act}} \cdot \text{Emb}\_{\text{Labor}} \tag{15}$$

where EmbFuel -MJ · l −1 and EmbLubr -MJ · l −1 are the embodied energy in the fuels and lubricants and EmbLabor (MJ · <sup>h</sup>−1) is the per hour labor energy. With respect to the indirect energy consumption, Equation (16) applies as follows:

$$\text{En}\_{\text{Indirect}} = \text{En}\_{\text{Manufacturing}} + \text{En}\_{\text{RM}} + \text{En}\_{\text{Houasing}} \tag{16}$$

The embodied energy for the manufacturing EnManufacturing (MJ) of the equipment (tractor and implement) is expressed in Equation (17), reduced to total hours of the operation:

$$\text{EnManfacturing} = \frac{\text{Emb}\_{\text{Manufacturing}} \cdot \text{M}}{\text{hife}} \cdot \text{t}\_{\text{act}} \tag{17}$$

where EmbManufacturing - MJ · kg−<sup>1</sup> is the embodied energy; M (kg) is the mass of the tractor and the implement, respectively and hlife (*h*) is the total life hours of the machinery. The repair and maintenance embodied energy EnRM (MJ), as already mentioned, is expressed as a percentage of the manufacturing energy requirements RMRate (%) (Equation (18)):

$$\text{En}\_{\text{RM}} = \text{En}\_{\text{Manufacturing}} \cdot \text{RM}\_{\text{Rate}} \tag{18}$$

Lastly, the housing energy EnHousing (MJ) is calculated by Equation (19) reduced to total hours of the operation as follows:

$$\text{En}\_{\text{Houseing}} = \frac{\text{Emb}\_{\text{Houseing}} \cdot \text{A}\_{\text{cover}}}{\text{h}\_{\text{life}}} \cdot \text{t}\_{\text{act}} \tag{19}$$

In that case, EmbHousing -MJ · <sup>m</sup>−<sup>2</sup> · <sup>y</sup>−<sup>1</sup> is the embodied energy of the facility required to house the equipment and Acover - m<sup>2</sup> is the area cover of each the machines (tractor or implement).

#### **3. Case Study Description**

The methodology described in the previous section was used to investigate the energy cost of performing two different tillage operations, considering different field areas and the available time window. The examined operations are secondary tillage performed by a field cultivator and deep tillage performed by a moldboard plow. Table 1 presents the operation parameters considered according to the ASABE standards [34]. The economic and total life is assumed to be identical for both implements. However, the operating speed differs with respect to the different operations because it depends on the operation depth, which is determined based on the general cultivation practices [43].


**Table 1.** Operation input parameters.

<sup>1</sup> [34], <sup>2</sup> [43] \* Commercial values.

As mentioned in the methodology, the estimated required implement width and the tractor power are adjusted to fit with the next bigger commercially available one. To that effect, feature values of commercially available machinery were collected for the calculations as they are provided by the manufactures. These specifications include the width, weight and dimensions of the implements and the power, weight, and dimensions of the tractors. All the tractors (with power ranging from 37.4 kW to 456 kW) that were inserted in the model are 4WD, while the fuel type is diesel in all the cases examined. Regarding the calculation of the transportation energy cost, it is assumed that each field is located at a distance of 1 km from the farm. The tractor travels at 40 km · <sup>h</sup>−<sup>1</sup> and performs each trip once in order to complete the operation and return to the farm.

The embodied energy parameters considered for the model calculation are presented in Table 2. Various methodologies have been proposed for the calculation of the manufactured embodied energy of farm machinery; however, for the purposes of this study, the approach of Kitani et al. was selected [35]. With respect to the embodied energy of the maintenance of the machinery, it is calculated as a percentage of the respective manufacturing energy according to the approach of Aguilera et al. [30]. The approach of Aguilera et al. was also used for the estimation of the housing embodied energy, considering the average value per year and per spatial coverage of the service building. All the parameters presented in this section represent the input parameters of the calculation process of Figure 2, which are highlighted in the green boxes.


**Table 2.** Energy parameters.

<sup>1</sup> [33], <sup>2</sup> [30], <sup>3</sup> [16], <sup>4</sup> [44].

#### **4. Results**

This section presents the results of the implementation of the energy calculation tool with respect to the estimation of the machinery energy cost of performing three different tillage operations. Two different assessment models were examined, namely the discrete and the continuous model. The discrete model refers to the calculation of the required energy using the adjustment to standard commercial implements and tractors according to the process presented in the methodology. The continuous model refers to the theoretical energy consumption as it would be estimated without the adjustment to standard implement and standard values. The machinery selected in the latter case would have the exact performance characteristics with the estimated requirements to perform the operation. Furthermore, in that case, there are no maximum available values, indicating that the operation is performed with one set of tractor–implement, irrespective of the width and power necessary to operate

within the time window available, which coincides with the actual operation time. The figures below present the total energy consumption, the total machine power, the number of machines, as well as the total operating width with reference to the time available and the cultivation area.

#### *4.1. Field Cultivator*

With respect to the secondary tillage with a field cultivator, Figure 3 presents the total energy cost, the total machine power, the number of machines, and the total operating width for cultivated areas of 100 ha in relation to the available time to perform the operation. The energy consumption for the specific area varies around the average value of 234.99 MJ · ha<sup>−</sup>1. The tool selects the required machinery system to perform the operation within the available time window. As expected, the number of machines required to complete the operation decreases as the available time increases. For example, three machines are required to complete the operation in 10 h and two machines when the available time exceeds 13 h. After 25 h, only one machine is required to complete the operation and, as a result, only the required width decreases. The total energy consumed shows an upward trend owing to increased total operation time until the 37 h available time. At this point, only the smallest implement is used to complete the operation, and the time of operation does not change further with the energy consumed, being the same irrespective of the time available.

**Figure 3.** Energy consumption, machine power, number of machines, and operating width for a cultivated area of 100 ha (field cultivator).

The discrete model represents the realistic approach for the calculation of the energy required to perform an operation. However, as presented in Figure 3, it is worth comparing it to the continuous model. In the case of field cultivator, the continuous model predicts an average energy cost of 111.83 MJ · ha−<sup>1</sup> for the area of 100 ha, which is 2.1 times smaller than the respective discrete. The tractor power and the implement width that are actually used are significantly larger than those that are ideally required, resulting in increased energy consumption.

Figures 4 and 5 present the energy cost and machinery requirements for completing the operation in 30 h and 50 h, respectively, in relation to the cultivated area. Little available time leads to a rapid increase of machinery needed, while in the case of more time availability, a second machine is required for fields larger than 200 ha. Considering that, for the field cultivator three, different implement options were available (4 m, 5 m, and 6 m), the operating width increases after an area of 135 ha for an available time of 50 h. The energy consumed shows a downwards trend when increasing the cultivated area, with average values of 218.94 MJ · ha−<sup>1</sup> for 30 h and 228.98 MJ · ha−<sup>1</sup> for 50 h available time. The average energy consumed increases when increasing the available time because the operating time increases as well. In the case of 50 h available, the energy consumed per hectare remains constant for fields between 2 and 135 ha with the given implements, as the smallest implement can cover this area without requiring a second one to perform the operation timely.

**Figure 4.** Energy consumption, machine power, number of machines, and operating width for 30 h time available (field cultivator).

**Figure 5.** Energy consumption, machine power, number of machines, and operating width for 50 h time available (field cultivator).

In order to gain full insight regarding the formulation of the model, Figure 6 presents the results for the time available and the area, ranging from 10 h to 100 h and from 10 ha to 300 ha, respectively. The discrete model estimates a higher energy consumption per ha than the continuous model. The only exception is the area of 1 ha, where the continuous model is more energy consuming than the discrete model after approximately 40 h of time available owing to the very small stipulated operating width. However, as such a small operating width is not realistic (approximately 0.015 m for 1 ha and 100 h available time), energy values corresponding to areas below 10 ha were discarded from Figure 6. The total average energy cost is 226.51 MJ · ha−<sup>1</sup> for the discrete model, which is 1.98 times higher than the continuous model. For the field cultivator, 69.4% of energy consumption derives from the operation of the machinery, while 30.6% is attributed to its ownership.

**Figure 6.** Energy consumption model for the field cultivator.

#### *4.2. Moldboard Plow*

Figure 7 presents the total energy consumption of performing primary tillage with a moldboard plow in a field of 100ha in relation to the time available to perform the operation. The average energy consumption for the 100 ha is 1505.82 MJ · ha−<sup>1</sup> for the discrete model. It is observed that the continuous model shows higher energy consumption, approximately 1.06 times, than the discrete model demonstrating the effect of the operation depth and type in the overall energy costs.

**Figure 7.** Energy consumption, machine power, number of machines, and operating width for a cultivated area of 100 ha (moldboard plow).

The required width and power to complete the operation on time do not differ significantly between the discrete and the continuous model. This can be attributed to the large energy requirements of the operation, as well as to the availability of implements. In the case of the moldboard plow, the commercially available width of the implements presented more subdivisions in the range of 1.1–4.4 m, providing the tool with more options during the machinery selection stage.

Figure 8 presents the evolution of the model with respect to the field area, considering time available of 50 h. For field areas up to 30 ha, the operation can be performed with one tractor carrying the minimum available implement (1.1 m). From this point, the implement width is further increased until it reaches the maximum available of 4.4 m. Then, with a further increase in the area, a second machine is needed to perform the operation timely. The width and power increasing pattern is the same for the second machinery as previously described until the point when a third machine is required. The average energy cost, in the case of 50 h of available time, is 1508.61 MJ · ha<sup>−</sup>1.

**Figure 8.** Energy consumption, machine power, number of machines, and operating width for 50 h time available (field cultivator).

The total energy consumption of the discrete and continuous model for moldboard plowing in a wide range of area and available time is presented in Figure 9. The continuous model has an average energy cost of 1604.95 MJ · ha<sup>−</sup>1, which is 1.06 times higher than the discrete model, which presents an average energy cost per hectare of 1510.37 MJ · ha<sup>−</sup>1. This finding highlights the importance of time when performing an agricultural operation of high operation energy requirements. For moldboard plowing, 90.4% of the energy consumed is attributed to operation, while only 9.6% of the energy consumed refers to the energy ownership costs. As a consequence, the adjustment to the closest bigger available implement resulted in a reduction of the total energy, as the actual time required to perform the operation decreased with the increased implements' width.

**Figure 9.** Energy consumption model for moldboard plow.

#### **5. Discussion**

In the previous section, the energy cost of performing secondary tillage with a field cultivator and primary tillage with a moldboard plow was presented. The results concerned tillage on medium soil with an operation depth of 15 cm and 30 cm, respectively. Tillage with a moldboard plow was calculated to be 6.67 times more energy consuming than field cultivation, which is an outcome that is expected considering the depth and the nature of the operation. For the field cultivator, the indirect energy costs account for almost 1/3 of the total energy consumed, while for the moldboard plow, only for 1/10. The direct energy cost is the determining factor of the total energy consumption. However, it is worth noting that its contribution to the total energy consumed increases as the energy intensity of the operation increases. Nevertheless, it should be stated that the available time to perform an operation and the different field areas do not seem to strongly affect the total energy per hectare.

Regarding machinery selection, the method employed focuses on the available time for field operations rather than the assessment of the financial cost of the equipment. In that light, the dimensioning of the implement requires sufficient power in the tractor as a prerequisite for the execution of the operation on time. The main disadvantage of the method is the fact that the operating economic costs are not considered during the selection stage. Additionally, the tool is designed in order to facilitate machinery selection, taking into account technical conditions rather than economic ones. As a consequence, a selection of the machinery system by applying this method can lead to non-optimum solutions in terms of economic criteria such as minimum operating or ownership cost. Furthermore, it is of limited use in cases where the equipment has already been purchased.

The results indicate that the proper allocation of tools and tractors and the use of the appropriate combination of machinery in each operation can affect the total energy consumed. In the case of light soil cultivation (field cultivator), the total energy consumed is higher than that expected from model calculations owing to the use of larger implements than those required to perform the operation on time. On the other hand, in the case of the moldboard plow, the use of a larger implement than the one required to perform the operation on time leads to a decreased energy consumption, which is justified when considering the contribution of the operation cost to the total energy cost of this operation. Nevertheless, the availability and the subdivision of the implement width strongly affect the diversion from the model consumption calculated for the examined cases.

Another factor that was examined was the sensitivity of the energy consumed with respect to the type of soil and the operation depth. In that end, the calculations were performed for different depths (plus 5 cm and minus 5 cm from the original example), as presented in Figure 10, and soil types (fine and coarse in addition to medium soil), as presented in Figure 11. For the field cultivator, it was found that the operation depth and the type of soil do not affect the total energy consumed in the case of the discrete model. This is attributed to the fact that the minimum available commercial machinery covers the requirements of the operation, irrespective of the type of soil or the depth.

**Figure 10.** Moldboard plow energy consumption for various soil types and operation depths.

However, in the case of moldboard plowing, both the operation depth and the type of soil strongly influence the total energy consumed. As presented in Figure 10, the energy consumed for deep

tillage in fine soil is almost two times higher than the that for coarse soil in the case of an operation depth of 25 cm. Additionally, in the case of fine soil, the energy for performing the operation at an operation depth of 35 cm is 2429.54 MJ · ha<sup>−</sup>1, 1.35 times higher than when performing it at a depth of 25 cm (1800.88 MJ · ha<sup>−</sup>1).

**Figure 11.** Moldboard plow energy consumption discrete model for various soil types.

#### **6. Conclusions**

The assessment and management of agricultural production is a multivariate problem. Thereby, the introduction and use of operation management tools can contribute towards the effective management of agricultural production and the quantification of its environmental impact. Existing tools and methodologies in most cases examine the entire life cycle of the cultivation or focus on the use of plant protection products and fertilizers owing to their proven severe adverse impact compared with the other agricultural inputs. Nevertheless, the quantification of the contribution of machinery to the total energy consumed is important and constitutes an essential first step for the evaluation of alternative technologies and tools. In fact, the results highlight the benefits of optimal allocation of tools, especially in the task of light soil cultivation, where the energy consumed using the commercially available machinery (226.51 MJ · ha−1) is almost double the optimal energy expressed by the continuous model (113.88 MJ · ha<sup>−</sup>1).

This paper attempted to examine and set the base of evaluating the energy consumption considering various operational parameters coupled with the commercial availability of implements and tractors. The model selects (from a collection of available implements and tractors) the required machinery to perform the operation in a given time, while evaluating the resulting energy cost. Even though the available time and the cultivated area do not seem to strongly affect the energy consumption per hectare, the commercially available implements highly affect the total energy cost. From the examined case studies, the effect of using unsuitable equipment was highlighted, indicating a promising potential for further reduction in the energy consumed from the machinery. Towards this direction, future work incorporating the consideration of financial parameters during the machinery selection is required in order to assess the economic and environmental performance of the machinery in an integrated manner.

**Author Contributions:** Conceptualization, M.L., D.K., C.G.S., and D.B.; data curation, D.K. and D.B; formal analysis, D.K.; investigation, M.P. and D.K.; methodology, M.L. and D.K.; software, M.L.; supervision, D.B.; validation, D.K., C.G.S., and D.B.; writing—original draft, M.L.; writing—review and editing, D.K., C.G.S., and D.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was supported by the project "Research Synergy to address major challenges in the nexus: energy-environment-agricultural production (Food, Water, Materials)" NEXUS, funded by the Greek Secretariat for Research and Technology (GSRT), Pr. No. MIS 5002496.

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

#### **References**


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

## **Decision Support System to Implement Units of Alternative Biowaste Treatment for Producing Bioenergy and Boosting Local Bioeconomy**

#### **Christos Vlachokostas 1,\*, Charisios Achillas 2, Ioannis Agnantiaris 1, Alexandra V. Michailidou 1, Christos Pallas 3, Eleni Feleki <sup>1</sup> and Nicolas Moussiopoulos <sup>1</sup>**


Received: 11 April 2020; Accepted: 30 April 2020; Published: 6 May 2020

**Abstract:** Lately, the model of circular economy has gained worldwide interest. Within its concept, waste is viewed as a beneficial resource that needs to be re-introduced in the supply chains, which also requires the use of raw materials, energy, and water to be minimized. Undeniably, a strong link exists between the bioeconomy, circular economy, bioproducts, and bioenergy. In this light, in order to promote a circular economy, a range of alternative options and technologies for biowaste exploitation are currently available. In this paper, we propose a generic methodological scheme for the development of small, medium, or large-scale units of alternative biowaste treatment, with an emphasis on the production of bioenergy and other bioproducts. With the use of multi-criteria decision analysis, the model simultaneously considers environmental, economic, and social criteria to support robust decision-making. In order to validate the methodology, the latter was demonstrated in a real-world case study for the development of a facility in the region of Serres, Greece. Based on the proposed methodological scheme, the optimal location of the facility was selected, based on its excellent assessment in criteria related to environmental performance, financial considerations, and local acceptance. Moreover, anaerobic digestion of agricultural residues, together with farming and livestock wastes, was recommended in order to produce bioenergy and bioproducts.

**Keywords:** bioenergy; efficiency of bio-resources; decision support system; multi-criteria analysis; sustainability

#### **1. Introduction**

Circular economy and resource efficiency have received increasing attention in research and environmental policy agenda in recent years. Circularity can be a catalyst for productive reconstruction and has a clear regional dimension, where the value of waste is a key element. Within the concept of circular economy, waste, energy, and water consumption need to be minimized [1–4] and, in any case, the corresponding environmental pressures should respect target or limit values [5]. Bioeconomy focuses on the use of biomass/biowaste in primary production procedures (e.g., agriculture, forestry, fisheries) and substantial valorization of raw materials. Biomass is expected to substantially support the achievement of the EU renewable energy targets [6].

Undoubtedly, there is a strong linkage between the bioeconomy, circular economy, bioproducts, and bioenergy. Streams of excess or end-of-life materials, previously regarded as waste, that originate from anthropogenic activities can be channeled through available technologies (e.g., anaerobic digestion, gasification, pyrolysis) and transmute into useable energy carriers, organic biofertilizer abundant in nutrients, and, in general, original and innovative materials. A characteristic example is the production of biofuel (biogas or syngas), which can be characterized as a major means of energy recovery from biowaste streams. An appreciable amplitude of different technological solutions is currently available for small, medium, or large-scale units of alternative biowaste treatment (UABT), with an emphasis on the production of bioenergy and other bioproducts. Such solutions are based on applications that extend from tailor-made systems, which incorporate existing available facilities (e.g., pumps and storages in an existing farm, buildings for the combined heat and power (CHP) installation etc.), to specialized concepts that are diverse, where the main parts have been pre-manufactured.

Although there are many solutions, biowaste disposal in landfills, or, in the worst case, in uncontrolled open dumps, is still an existing practice internationally [7]. This cannot be considered as a proper managerial approach for biowaste, having in mind the intense environmental load imposed. Apart from aesthetic degradation, impacts comprise air pollution and Green House Gas (GHG) emissions, soil and water contamination, reduced land values, and landscape blight. Consequently, optimal biowaste management constitutes a critical activity that reinforces environmental sustainability by minimizing impacts, especially related to uncontrolled discarding [8–12].

The appropriateness and benefits of each technological solution should be examined carefully to support treatment decisions, considering the type of biowaste, local characteristics, and conditions [13–15]. In the material to follow, a generic methodological scheme was developed and demonstrated to support optimal decision-making in promoting and implementing UABTs. The basic structure and components combine environmental, economic, and social parameters and multi-criteria decision analysis (MCDA) for the optimal site selection of a UABT. Since the early 1980s, Ross and Soland [16] reported that problems involving the location and promotion of such units include multi-criteria considerations and ought to be modelled likewise. Nevertheless, the literature has scarcely studied special waste streams, such as biowaste, centralizing mostly on municipal solid waste (MSW) management [17–22]. MCDA has gained wide acceptance in recent years over quantitative modeling, as MCDA embodies both quantitative as well as qualitative variables [23] and can be tractably combined with other tools, such as life cycle assessment, ecological footprint, and environmental indicators [24]. The approach includes the social criterion (e.g., NIMBY - Not In My Back Yard syndrome), which is usually disregarded [25,26]. It also assists relevant stakeholders, public authorities, and producers by providing a roadmap to the feasibility and essential steps for the development of a micro-to-medium-scale localized UABT. Localization is a prerequisite to achieve economies of scale and an emphasis is given to an optimal location, and the production of bioenergy and other bioproducts [27,28].

In this context, the present work frames a consistent generic methodology and seeks optimal UABT solutions, based on multiple criteria and parameters. More specifically, apart from defining the optimal site for the installation of a new UABT, the proposed methodology also assesses the recommended technological solution, based on the characteristics and quantities of locally available biowaste, as well as its size and capacity. The methodology is demonstrated for a real-world case in the region of Serres, Greece.

#### **2. Methodology: Basic Structure and Components**

#### *2.1. Selection of Area and Inventory Analysis*

Figure 1 presents the main components and the basic steps of the methodological scheme. More specifically, the first step is to determine the wider area for the development of a UABT (Step 1). Crucial aspects for this decision are the availability, quantity, and quality of biowaste regionally. Characteristic typologies can be found in the Waste Framework Directive [29] and related regulations, e.g., park and yard waste, animal waste (feces, urine, and manure), kitchen and food waste, vegetal waste (food preparation and products), household and similar wastes, common sludges, mixed waste, undifferentiated materials, etc. To assure smooth replenishment and long-term sustainability of a UABT, it is important to pre-assess the continuing availability and seasonality of biowaste locally. Regional logistics infrastructure, mainly transportation and storage facilities, also need to be considered for the available feedstock, due to their criticality in the unit's viability, since the transportation cost is a crucial parameter for the viability of the logistics network [30]. All the above data can be the output of a strengths, weaknesses, opportunities, and threats analysis, for the area under consideration. This initial stage of the methodology puts forward a pool of available types of biomass as potential feedstock for an UABT.

**Figure 1.** Basic structure and components of the methodological framework.

In the case that an area is considered "eligible" for the development of a UABT, the second step is the preparation of a detailed biowaste inventory analysis that (i) identifies meticulously the biowaste sources, (ii) calculates the available biowaste quantities, and (iii) describes the characteristics of the available biowaste (Step 2). Such an inventory will detail the availability of exploitable biowaste within the region under study and will assist decision-making.

#### *2.2. Optimal Selection of the Site*

The inventory analysis is followed by the multiple criteria decision analysis (MCDA) for the optimal location of the UABT (Step 3a). This process is strongly interrelated and usually realized in parallel to the selection of the most appropriate technology (Step 3b). Multi-criteria mathematical modeling examines criteria that are usually in conflict in the decision-making process [31–36]. On this basis, and after screening all legislative constraints (e.g., environmental permits), the appropriate multi-criteria analysis technique needs to be determined. The literature reports a considerable number of techniques available, with different characteristics and uses [37,38]. Multi-criteria analysis techniques are more or less suitable depending on the special characteristics of the case under consideration [39]. In the scientific literature, modelling of waste management solutions mostly takes into account the economic and environmental dimensions of the problem, whilst the social concerns are usually neglected. Morrissey and Browne [26] meticulously reported sustainable waste management applications. It should be noted that the social pillar constitutes a crucial matter for the future viability of the investment in cases like the one examined in the material to follow. On this basis, social concerns should be simultaneously considered for deciding on the optimal UABT site.

The methodology developed herein adopts the ELECTRE III technique [40]. A main characteristic of this technique is the use of the pseudo-criteria concept in order to depict the different angles of the studied case/problem. ELECTRE III utilizes three pseudo-criteria and initiates a robust comparison of each option with another in relation to all criteria included in the analysis. An ascending and descending distillation process forms the basis for the construction of two complete pre-orders. The intersection of the two complete pre-orders produces a final classification of the alternatives. An advantage of this technique is the sensitive analysis capability. In the framework of the sensitivity analysis, scenarios for the values of the main parameters are performed in order to further test the solution by observing the effect of the adopted values' variation to the result. Comparative assessment of the pre-orders for each scenario leads to a final robust outcome or to a model re-analysis [41]. This technique is widely used in the examination of environmental problems and waste management issues, which provides an advantage to the others [20,33,42–44].

The designation of three thresholds is a prerequisite, i.e., (i) preference threshold (p), (ii) indifference threshold (q), and (iii) veto threshold (v). ELECTRE III uses these thresholds in its mathematical rationale in order to better include real-life uncertainties [41]. Another advantage of this multi-criteria technique is the simultaneous consideration of quantitative (e.g., price of sites, distances, etc.) and qualitative criteria (e.g., landscape degradation, aesthetics, social acceptance, etc.), as this approach depicts a good fit of the data in such applications.

It goes without saying that a meticulous investigation in the wider area under study is needed in order to landmark all the available locations appropriate for the installation of a UABT. As a next step, the assessment of different, usually conflicting, criteria is performed in parallel with defining the weighting factor of each criterion. The weighting factor expresses its relative significance in comparison to the others. A clear definition of the parameters included in the analysis is necessary in order to reliably value all available alternatives and perform the pairwise comparison process. After defining the criteria and weighting factors, the related information and data should be assembled. Multi-criteria evaluation of sites for the UABT location is mathematically formulated with the use of a set of criteria (Cr1, Cr2, Cr3 ... ) applicable to a set of alternatives (A1, A2, A3 ... ). Vj(A) mathematically expresses the evaluation of alternative A i, for the criterion j. ELECTRE III is a ranking method, thus it puts forward a ranking prioritization, which is grounded on binary outranking relations for two interrelated concepts, i.e., (i) concordance (cj) and (ii) non-discordance (dj). More specifically, the concept of the concordance relation is applicable when alternative A1 outranks alternative A2 in case a sufficient majority of criteria are in favor of alternative A1. The non-discordance relation is applicable when the concordance condition holds, and none of the criteria in the minority should be opposed strongly to the outranking of A2 by A1. The credibility index characterizes the assertion that A1 outranks A2 and illustrates the true degree of this assertion [45]. A pair of alternatives is compared for each criterion with the use of pseudo-criteria, namely the indifference (qj) and preference (pj) for which they apply:


The concordance index cj (A1, A2) of each criterion j is mathematically formulated as follows:

$$\begin{aligned} \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{1}) - \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{2}) &\leq \mathbf{q}\_{\mathbf{j}} \Leftrightarrow \mathbf{c}\_{\mathbf{j}}(\mathbf{A}\_{1}, \mathbf{A}\_{2}) = 0 \\ \mathbf{q}\_{\mathbf{j}} &< \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{1}) - \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{2}) < \mathbf{p}\_{\mathbf{j}} \Leftrightarrow \mathbf{c}\_{\mathbf{j}}(\mathbf{A}\_{1}, \mathbf{A}\_{2}) = \frac{\mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{1}) - \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{2}) - \mathbf{q}\_{\mathbf{j}}}{\mathbf{p}\_{\mathbf{j}} - \mathbf{q}\_{\mathbf{j}}} \\ \mathbf{V}(\mathbf{A}\_{1}) - \mathbf{V}(\mathbf{A}\_{2}) &\geq \mathbf{p}\_{\mathbf{j}} \Leftrightarrow \mathbf{c}\_{\mathbf{j}}(\mathbf{A}\_{1}, \mathbf{A}\_{2}) = 1 \end{aligned}$$

A global concordance index CA1A2 for each pair (A1, A2) is calculated with the use of cj (A1, A2) as mathematically illustrated below:

$$\mathbf{C}\_{\mathbf{A}\_1\mathbf{A}\_2} = \frac{\sum\_{\mathbf{j}=1}^{\mathbf{n}} \mathbf{w}\_{\mathbf{j}} \cdot \mathbf{c}\_{\mathbf{j}}(\mathbf{A}\_1, \mathbf{A}\_2)}{\sum\_{\mathbf{j}=1}^{\mathbf{n}} \mathbf{w}\_{\mathbf{j}}},\\\mathbf{w}\_{\mathbf{j}} \text{ expresses the weighting factor of criterion j.}$$

The discordance index (dj) is computed with the adoption of indifference (qj), preference (pj), and the veto threshold (vj), which expresses the maximum acceptable difference for not rejecting the assertion A1 outranks A2. More specifically:


Thus, dj (A1,A2) is mathematically formulated as:

$$\begin{cases} \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{2}) - \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{1}) \leq \mathbf{p}\_{\mathbf{j}} \Leftrightarrow \mathbf{d}\_{\mathbf{j}}(\mathbf{A}\_{1}, \mathbf{A}\_{2}) = 0 \\\ \mathbf{p}\_{\mathbf{j}} < \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{2}) - \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{1}) < \mathbf{v}\_{\mathbf{j}} \Leftrightarrow \mathbf{d}\_{\mathbf{j}}(\mathbf{A}\_{1}, \mathbf{A}\_{2}) = \frac{\mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{2}) - \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{1}) - \mathbf{p}\_{\mathbf{j}}}{\mathbf{v}\_{\mathbf{j}} - \mathbf{p}\_{\mathbf{j}}} \\\ \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{2}) - \mathbf{V}\_{\mathbf{j}}(\mathbf{A}\_{1}) \geq \mathbf{v}\_{\mathbf{j}} \Leftrightarrow \mathbf{d}\_{\mathbf{j}}(\mathbf{A}\_{1}, \mathbf{A}\_{2}) = 1 \end{cases}$$

.

The index of credibility δA1A2 of the claim A1 outranks A2 is then formulated as:

$$\delta\_{\mathbf{A}\_1 \mathbf{A}\_2} = \mathbb{C}\_{\mathbf{A}\_1 \mathbf{A}\_2} \prod\_{\mathbf{j} \in \bar{\mathcal{F}}} \frac{1 - \mathbf{d}\_{\bar{\mathbf{j}}}(\mathbf{A}\_1, \mathbf{A}\_2)}{1 - \mathbb{C}\_{\mathbf{A}\_1 \mathbf{A}\_2}}, \text{ with } \bar{\mathcal{F}} = \langle \mathbf{j} \in \mathcal{F}, \mathbf{d}\_{\bar{\mathbf{j}}}(\mathbf{A}\_1, \mathbf{A}\_2) > \mathbb{C}\_{\mathbf{A}\_1 \mathbf{A}\_2} \rangle.$$

The claim A1 outranks A2 is rejected if vj is exceeded for at least one of the criteria under consideration. After constructing the ranking scheme, as a last step, sensitivity analysis is activated, considering that values and assessments are often subjective in real-world cases and originate from less or more reliable estimations (criteria qualitative expression, threshold parameters, weighting factors, etc.). As already highlighted, this is considered a major advantage of the adopted MCDA technique due to the high data uncertainty in the thematic area under study [46].

#### *2.3. Technology, Size, and Capacity of the UABT*

The decision for the optimal site in the wider area presents a strong interrelationship with the determination of the technology to be adopted (Step 3b). A range of available technologies are available (e.g., anaerobic digestion, aerobic digestion, gasification, pyrolysis) and possibly a combination of these. The key factors that need to be considered for selecting the optimal technology are the potential for bioenergy production (electricity, heat or gas), possible utilization of the residues of each technological solution (e.g., biofertilizer), suitability of available bioresources, and pre-treatment requirements. All this information supports the optimal decision on technology for the UABT, based on constraints related to the type of biowaste and specific technological solutions used.

As a next step, the assessment of the size and capacity of the UABT is realized (Step 4), which is based on a cost-benefit analysis (CBA) and the determination of the biomass mixture for the UABT (range of sources depending on the biowaste quality). The latter ultimately leads to the dimensioning and detail of the technical specifications for the UABT. Highlighting strategies for risk management is the last step of the methodological roadmap presented (Step 5). The abovementioned methodological framework is expected to support policy-makers to facilitate the conceptual phase towards the development of a UABT, before advancing to the next phase of the planning phase where more resources are required. In this light, decision-makers may screen available options and give a green light for further analysis only to those cases where the endeavor is proven viable.

#### **3. Application of Methodology in the Region of Serres, Greece**

#### *3.1. The Area under Study*

The applicability of the methodological framework was validated in a real-life case study in the region of Serres, Greece. The region of Serres under study was pre-selected based on the considerable regional agricultural activity and biomass availability. The population of the wider Serres area is approximately 200,000, out of which 25% is urban, 20% semi-urban, and 55% rural/agricultural. Its economic activity mainly deals with the primary sector and especially with agriculture (approximately 60%). In the region, there are 162,800 ha of cultivated land and 113,300 ha of pasture land. More than 90% of the cultivated land is private. In the area, there is well-developed agricultural production, which contributes 3.4% of the total national agricultural product, with the main drivers being dairy and meat production. Apart from farmers, within a range of 20 km from the municipality of Serres, there are a number of agri-food manufacturers. Specifically, the relevant production activity consists of 2 leading dairy industries, 6 small cheese/yogurt production units, 2 olive mills, 1 potato packing and processing unit, 4 poultry farms, more than 30 cattle and pig rearing units for dairy and meat production, and one slaughterhouse.

These facilities are capable of supplying adequate quantities of biowaste to generate a mixture that ensures the continuous operation of a UABT. Considering the available data and given the population and basic economic activity characteristics of the wider area under study, a primary course of action regarding biowaste energy content could be supported by the recent legislative framework set in Greece regarding energy communities [47]. The legislation sets the rules for the incorporation of quasi-private entities with the option of participation of the municipalities as shareholders, while also providing the alternative of subsidized (or with tax-incentives) renewable energy plants. Such strategic schemes could also opt from the formation of public–private partnerships.

The region under study has great potential for the application of a productive UABT model due to many reasons, such as: (i) The availability of resources and biowaste, (ii) scientific potential and know-how, and (iii) the existence of a primary sector with growth potential and needs for modernization and reduction of production costs. Last, but not least, local farmers are strongly interested in adjusting their traditional economic attitudes. This high level of engagement was is result of intense awareness-raising activities (including flyers/brochures dissemination, press and media publicity, public presentations, visits to demonstration projects) as well as public opinion surveys that leaded to accurate expectation mapping of all involved local stakeholders [48].

#### *3.2. Alternative Sites and MCDA*

According to the methodology described, different locations for the development of the UABT within the area under study (Figure 1) were promoted. Following a detailed survey on land availability, four alternative available fields were selected, such as those sites are summarized in Table 1.

Next, the criteria for the analysis were defined in order to interlace financial viability, environmental performance, and acceptance by the local community, taking additionally into account all legislative constraints. A critical mass of local experts, also aware of the local special characteristics, were invited. In total, 13 experts representing local and regional authorities and the business sector in the study area were interviewed to decide on the criteria to be used in the case under study. This resulted in the following criteria: (Cr1) Distance from biowaste source, (Cr2) accessibility of the site with the existing infrastructure, (Cr3) distance from units for the consumption of energy (electricity and heat) produced, (Cr4) distance from farms for compost use, (Cr5) value of land, (Cr6) impacts on local populations (e.g., air pollution, odors, noise pollution), (Cr7) impacts on local ecosystems, (Cr8) aesthetic degradation, (Cr9) increase in road traffic due to UATB's operation, and (Cr10) social acceptance.


**Table 1.** Alternative locations for UABT development in the area under study.

The four site locations were assessed over their performances on the 10 designated criteria (Table 2). Criteria Cr2, Cr6, Cr7, Cr8, Cr9, and Cr10 were qualitatively assessed by the experts involved in the survey. Those criteria that are related to the distances (Cr1, Cr3, Cr4) were assessed with the use of Google maps and the location of available biowaste sources and facilities where the outputs of the UABT could be exploited (farms and manufacturers) [48]. Cr5, which is related to the value of the land, was assessed with mean objective values, which were provided by the Hellenic Statistical Authority. For all selected criteria, the higher the performance, the more preferable the alternative is assessed.


**Table 2.** Assessment of alternative site locations (*in 1–10 scale, where 10 is the most preferable*).

Weighting factors are also depicted in Table 2. Those were estimated as averages of the experts' views. In order to overcome subjectivity issues, sensitivity analysis was employed. In our case, a low computational time was required to re-calculate the optimal solutions with the modified parameters. The calculation of pi (preference thresholds) for the selected criteria was based on the use of Equation (1) [44,46,49,50], while the calculation of qi (indifference thresholds) was based on the use of Equation (2) [46,51]:

$$\mathbf{p}\_{\mathrm{i}} = \frac{1}{\mathbf{n}} (\mathbf{V\_{i}} \mathbf{m} \mathbf{a} \mathbf{x} - \mathbf{V\_{i}} \mathbf{m} \mathbf{n}),\tag{1}$$

$$\mathbf{q}\_{\rm i} = 0.3 \cdot \mathbf{p}\_{\rm i} \tag{2}$$

#### *3.3. Technology, Size, and Capacity of the UABT*

Based on the above analysis, a combination of anaerobic and aerobic digestion is promoted as the selected technological solution. Anaerobic digestion (AD) fits the local agricultural and farming well since manure, cheese whey, side products (i.e., rotten potato pulp, oil mill waste), and slaughterhouse waste are efficient substrates available in the area. AD is considered as the process where value is produced mainly for the generation of bioenergy (CHP) and biofertilizer. AD is a fermentation process, which is realized in an air-sealed biodigester (under conditions of oxygen absence). The feedstock (biowaste substrates) is converted into fuel gas and digestate (as a bio-product). The biogas that is produced comprises mainly of CH4 (50%–70%) and CO2 (30%–50%). Moreover, in smaller quantities, the biogas consists of H2O vapor, H2S, and other elements. The digestate is produced from the digestion of substrates, after biogas extraction. Wet digestate comprises of nutrients (e.g., nitrogen, lignin, phosphorous), inorganic salts (e.g., ammonium, phosphate, potassium), as well as other minerals. Therefore, it is suitable for use as a naturally derived fertilizer. In the literature, there are many studies related to the exploitation of AD [52–55]. Regional AD applications are also detailed in [56]. The substrates used determine the appropriate technical infrastructure that is required (e.g., pipes' sizing, pumps' specifications, gas storage requirements, gas treatment technology, and Combined Heat and Power - CHP design).

The technology adopted does not add to the CO2 load in the atmosphere, since the CO2 produced is offset by the prevented emissions of CH4 in the case that slurry is stored in open spaces. In this light, biogas can significantly contribute to decarbonization of the economy. Furthermore, decentralized bioenergy generation from biofuel in UABTs with CHP is becoming more attractive for the cases where generated heat is exploited in facilities within the proximity of the UABT. Moreover, the use of UABTs' end-products (energy and bioproducts) needs to be meticulously examined within the pre-planning phase. In this context, the use of the produced bio-fertilizer in farms close to the developed UABTs highly supports their financial performance and effectiveness. The sufficiency of arable land for spreading available biofertilizer and an available adjacent market for the digestate needs to be investigated. In addition, biofertilizers' nutrients are predominately contained in a mineral form, which allows ease of absorption by plants and increased uptake efficiency, in comparison to nutrients in raw manure or in slurry [57].

Based on the above, it is proven in practice by the operation of such multi-functional facilities that this technology offers reduced local community energy costs, low-cost and environmentally safe recycling of manure and biowaste, cheap and efficient crop biofertilization means, and a reduction of odors due to intensive farming activity. Thus, considering also that in many areas, organic substances are still disposed in landfills (contributing to local CH4 emissions), the proposed UABT apparently dually contributes towards global warming mitigation through (a) the decreased CH4 emissions in animal farms (due to digestate management and avoidance of open slurry storage), and (b) the production of a green decarbonized (bio-)fuel. Along with the latter, improved digestate nutrient management, in combination with sustainable agri-practices [58], has a positive effect on the reduction of NH3 and NOx emissions, as well as the eutrophication of surface and ground water since leakages can be decreased or even avoided. Furthermore, local eutrophication is also achieved through biogas treatment of manure [59].

#### **4. Results and Discussion**

#### *4.1. Feedstock*

Based on the inventory analysis, the UABT receives organic biowaste types. More specifically, Table 3 presents the results of the annual received quantity of each substrate, their supply (d/a), and the quantity of each substrate per day of supply [48]. The codification of each substrate is presented in accordance with the European Waste List (EWL) (as per Annex to Decision 2000/532/EC).


**Table 3.** Substrate mix composition, yearly supply schedule, and yearly supply quantities.

The energy density of manure is considerably low, mainly due to the high content in water, along with the relatively low gas yield. The latter decreases its attractiveness for the case of long transportation distances. On the other hand, hydraulically, slurry is comparatively easy to handle. In the cases where the proportion of manure is relatively high within the digester, substrates like grass or solid manure, which are usually hydraulically more demanding, show a good performance in a small-sized biogas plant.

#### *4.2. Optimal Site*

The use of ELECTRE III in the case herein examined was realized with LAMSADE software. The two distillations (ascending and descending) that result from the analysis of the sites' performances are graphically displayed in Figure 2a. Both ascending and descending distillations promote Site 2 as the optimal alternative for the UABT's location. This is also displayed in Figure 2b, where the four alternatives ae hierarchically ranked. Site 2 (optimal location) is approximately 1.5 km south-southwest of the city of Serres, within municipality-owned land of approximately 16 acres.

As presented in Figure 2b, site 2 is judged as the optimal location locally. The latter is grounded on its outstanding assessment in most of the examined criteria (Table 2). This is well-justified due to the fact that site 2 is in the proximity of major producers of high-value biomass and potential end-users of the produced energy. Due to the fact that the area is industrial, the infrastructure in the proximity of site 2 is in very good condition, while the specific site also shows a very good performance in the social criteria (impacts on local populations and ecosystems, aesthetic degradation, increase in road traffic, social acceptance).

As a next step, a sensitivity analysis was conducted. The thresholds in the mathematical formulation (preference and indifference) were modified, as presented in Table 4. More specifically, six scenarios were examined. The optimal location for the development of the UABT in the region of Serres is site 2 for all examined scenarios. The latter increases the robustness of the basic solution and provides confidence in relation to the efficiency of the location over other available alternatives.

**Figure 2.** Optimal UABT location: (**a**) Ascending and descending distillations - (**b**) Ranking of alternative locations.


**Table 4.** Sensitivity analysis.

#### *4.3. Technical Specifications*

The biowaste is collected from the production and disposal sites and transported by container trucks and/or simple trucks to the UABT. Biowaste is introduced into the reception tanks to be mixed through agitation and formulate the final substrates' mix for AD. The mix is pumped and transferred (via the displacement pumps) to the digesters. The UABT will produce approximately 160,000 m<sup>3</sup> of biofuel per year. The biogas is collected and transferred to the internal combustion engine (ICE) for the production of electrical and thermal energy. The estimated electricity generation amounts to 311,996 kWh/a of operation, including self-consumption needs and grid losses. The estimated thermal bioenergy generation is 356,216 kWh/a, including thermal losses and self-thermal energy consumption. More specifically, the excess electrical energy is fed to the low voltage grid and its selling price offsets the electricity purchase costs of the municipality of Serres. The excess thermal energy serves the heating needs of a greenhouse adjacent to the unit's site and is partially used to heat the substrates' mix inside the biodigesters through a shell-type heat exchanger. This conducts the heat from the ICE's exhaust gases to the closed hot water circuit passing through the digesters. At the closed hot water circuit's return, the water, at a lower temperature than that when coming out of the shell-type heat exchanger, is pre-heated with the use of a heat exchanger (plate-type) and simultaneously cools the ICE, before passing through the shell-type heat exchanger again and increasing its temperature from the exhaust gases of the ICE.

The above process was calculated to produce 3121 t of aqueous residue, rich in nutrients for deposition in cultivations adjacent to the UABT. The proposed AD process is thermophilic. The biowaste mixture remains inside the digesters for at least 20 days (hydraulic retention time ≥ 20 days) at a temperature of at least 52 ◦C. The digestate is pumped in a controlled manner and transported to the screw separator, which recovers its solid fraction to be bagged, while the liquid fraction is placed in a lagoon. The solid fraction bags are picked up by trucks and the liquid fraction is pumped and transported by tank containers, both to be deposited in adjacent crops during eligible fertilization periods.

It should be emphasized that the application of biofertilizer needs to be realized in the growth season for any given crop. This is necessary in order to increase the uptake of nutrients, as well as to avoid a leaching or surplus of nutrients. Apart from the nutrient content, a fertilizer needs to also show high quality in sanitation and the reduction of pathogens; quantities of organic compounds and heavy metals, plastic, or other contaminants; etc.

The proposed methodology's clear strength lies on the multi-functionality of this technological solution. UABTs' sustainability can be grounded on many different aspects, such as (a) exploitation of organic waste, (b) upgraded management of manure, (c) increased uptake efficiency of nutrients, (d) reduced odors, (e) protection of the environment and reduced GHG emissions, (f) increase of material value, and (g) bioenergy (in the form of electricity and heat) and biofuel production.

#### **5. Conclusions**

In our work, we proposed a methodological scheme aimed towards supporting decision-makers in assessing the development of a UABT, based also on multiple criteria and factors. The methodology was effectively validated in a real-world case study in the region of Serres, Greece. In the approach proposed, the decision for the optimal location of the UABT within the area under study was based on the results of a multicriteria analysis (with the use of ELECTRE III technique), taking into consideration qualitative and quantitative criteria. In the case presented, the optimal site for the location of the UABT results from its excellent performance in the selected criteria, which represent all pillars of sustainable development, namely financial viability, public acceptability, and environmental protection.

The proposed methodological scheme provides a roadmap either for public bodies or private companies that aim to invest in green energy projects. The methodology followed in this work is generic so as to be followed under different requirements and constraints in any real such investment, with slight modifications of the criteria or values in the thresholds and weighting factors according to the special conditions in the area under study. The methodological framework can also be used to make decisions on other issues related to the supply chain management, namely the optimal location of collection sites, warehousing, sorting centers, etc., after relevant modifications in the criteria considered.

Even though the production of biogas has been widely realized over the past decades, the large-scale substitution of fossil fuels is still considered a major innovation in the energy sector. Energy technologies that exploit biomass for energy production are still non-competitive and technologically inefficient and immature against the conventional use of non-renewable sources. The methodology presented provides a little stepping stone towards supporting a bioeconomy within the field of clean technology biowaste and residues, within the concept of a circular economy, which is currently in its advent. It should be emphasized that such AD technological solutions are followed by challenges that limit the wider applicability, i.e., (i) the lack of markets for the digestate produced resulting from poor public perception; (ii) lack of recognition of AD and composting, especially regarding nutrient recirculation; (iii) purity of segregated biowaste streams; and (iv) requirements in capital investment.

**Author Contributions:** Conceptualization, methodology and paper writing, C.V.; Conceptualization, graphical illustration and validation of results C.A.; AD-related data, I.A.; Software runs and sensitivity analysis, A.V.M.; Case study related data, C.P.; Literature review, E.F.; Structure formulation, N.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the INTERREG IPA Cross Border Cooperation Programme CCI 2014 TC 16 I5CB 009 grant number SC032 (ZEFFIROS project).

**Acknowledgments:** The authors would like to wholeheartedly thank the experts that participated to the MCDA application in the study area.

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

#### **References**


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

#### *Article*

## **Forecasting of Day-Ahead Natural Gas Consumption Demand in Greece Using Adaptive Neuro-Fuzzy Inference System**

#### **Konstantinos Papageorgiou 1,\*, Elpiniki I. Papageorgiou 2,3, Katarzyna Poczeta 4, Dionysis Bochtis 3,\* and George Stamoulis <sup>5</sup>**


Received: 31 March 2020; Accepted: 4 May 2020; Published: 7 May 2020

**Abstract:** (1) Background: Forecasting of energy consumption demand is a crucial task linked directly with the economy of every country all over the world. Accurate natural gas consumption forecasting allows policy makers to formulate natural gas supply planning and apply the right strategic policies in this direction. In order to develop a real accurate natural gas (NG) prediction model for Greece, we examine the application of neuro-fuzzy models, which have recently shown significant contribution in the energy domain. (2) Methods: The adaptive neuro-fuzzy inference system (ANFIS) is a flexible and easy to use modeling method in the area of soft computing, integrating both neural networks and fuzzy logic principles. The present study aims to develop a proper ANFIS architecture for time series modeling and prediction of day-ahead natural gas demand. (3) Results: An efficient and fast ANFIS architecture is built based on neuro-fuzzy exploration performance for energy demand prediction using historical data of natural gas consumption, achieving a high prediction accuracy. The best performing ANFIS method is also compared with other well-known artificial neural networks (ANNs), soft computing methods such as fuzzy cognitive map (FCM) and their hybrid combination architectures for natural gas prediction, reported in the literature, to further assess its prediction performance. The conducted analysis reveals that the mean absolute percentage error (MAPE) of the proposed ANFIS architecture results is less than 20% in almost all the examined Greek cities, outperforming ANNs, FCMs and their hybrid combination; and (4) Conclusions: The produced results reveal an improved prediction efficacy of the proposed ANFIS-based approach for the examined natural gas case study in Greece, thus providing a fast and efficient tool for utterly accurate predictions of future short-term natural gas demand.

**Keywords:** neuro-fuzzy; ANFIS; neural networks; soft computing; fuzzy cognitive maps; energy forecasting; natural gas; prediction

#### **1. Introduction**

The increasing technological advancements and the rapid global population growth have led to a remarkable increase in energy consumption all over the world and especially in the developed and developing countries. After all, energy consumption is an index of a society's economical welfare and represents the economic development of a city or country [1]. Due to this unexpected boost in energy consumption over the past few decades, the need for energy demand management became crucial for achieving economic success that will result in self-sufficiency and economic development [2]. Thus, energy consumption forecasting is essential, as it predicates an energy efficient policy, the optimization of usage [3] and energy supplies optimum management. Even though a variety of methods have been investigated for energy demand forecasting, this is not an easy task, as it is affected by uncertain exogenous factors such as weather, technological development and government policies [4].

Among all energy resources, natural gas (NG) has received the largest increase in consumption lately [5], mostly due to its popularity as a clean energy source, with respect to environmental concerns. In particular, this energy source is characterized by low-level emissions of greenhouse gases in comparison with other non-renewable energy sources [6,7] and is considered as the cleanest-burning fossil fuel [8–11]. One fifth of the world's primary energy demand is covered by NG [7] and is linked to industrial production, transportation, health, agricultural output and household use.

According to the sharp increase in NG consumption, the forecasting of NG consumption has attracted great attention since it is essential for project planning, gas imports, tariff design, optimal scheduling of the NG supply system [12], indigenous production, infrastructures planning and cost reduction at different levels [7]. Moreover, some noteworthy factors like the need for distribution planning, especially in residential areas, the increasing demand of NG and the restricted NG network in many countries, make consumption forecasting on an hourly, daily, weekly, monthly or yearly basis highly important [13].

Especially in national energy strategy, NG demand forecasting is of high importance and can help policymakers all over the world to choose certain strategies in this direction. The development of systems that model NG consumption could assist in good government policymaking. Over the past decade, there has also been a significant increase in NG consumption in Greece and so the prediction of demand has become crucial accordingly. In particular, there is a need for further distribution planning in the Greek territory, especially in residential areas and in cases of high demand, when the accumulation ability of the network itself is decreased [13].

Selecting appropriate load forecasting methods is the most important step. Various models have been proposed in the field of energy forecasting so far, and can be divided into three groups: (i) traditional statistical models such as regression analysis, time series methods and ARIMA [14]; (ii) artificial intelligence (AI)-based methods such as wavelet analysis, artificial neural network (ANN) methods, neuro-fuzzy, machine learning, gray theory prediction; and (iii) hybrid models, which are a combination of various forecasting methods [15–23]. As regards the forecasting period horizon, prediction of demand is classified into three categories: short-term, medium-term and long-term forecasting. Short-term forecasting is used for hourly, daily and weekly demand predictions and in regard to NG, it is oriented mainly in system management and balancing, storage capacities optimization, as well as in making optimum purchasing and operating decisions [24,25]. Medium-term forecasting deals with seasonal demand prediction (one to several months) and mainly plans the fuel purchases, while long-term forecasts (more than a year ahead) aim to develop the power supply and delivery system [26,27].

ANN, genetic algorithm (GA) and fuzzy inference systems (FIS), as artificial intelligence (AI) techniques, are among those methods that are often used in the energy domain and specifically for energy demand forecasting, due to their high flexibility and reasonable estimation and prediction ability. From the relevant literature, there have been several attempts in energy forecasting by ANNs, like those in [28–32]. ANNs have been applied to forecast electric energy consumption in Saudi Arabia [33], energy consumption of a passive solar building [34], energy consumption of the Canadian residential sector [35], the peak load of Taiwan [36], while ANNs have been further explored in [37–41] for short-term load forecasting. An abductive network machine learning for predicting monthly electric energy consumption in domestic sector of Eastern Saudi Arabia was proposed in [42]. Moreover, support vector machines (SVMs) and genetic algorithms (GA) were explored in [29,30,43] to predict

electricity load. In [44], SVMs coupled with empirical mode decomposition were used to perform long term load forecasting. GAs have been used in [45] and [46] to estimate Turkey's energy demand and electricity demand in the industrial sector, respectively. In addition, Azadeh et al. have proposed the integration of GAs and ANNs to estimate and predict electrical energy consumption [47]. A number of literature reviews have also been performed regarding various energy fields. For example, in [48], a review of the conventional methods and AI methods for electricity consumption forecasting was provided, while in [49] the strengths, shortcomings, and purpose of numerous AI-based approaches in the energy consumption forecasting of urban and rural-level buildings were discussed. In [50], a review about conventional models, including time series models, regression models and gray models was conducted with respect to energy consumption forecasting. Moreover, an overview of AI methods in short term electric load forecasting area was discussed in [51].

#### *1.1. Indicative Related Work on AI Applied in Natural Gas Consumption Forecasting*

On the other hand, there are many studies where different AI methods like neural networks, neuro-fuzzy and other ANN topologies have been investigated and applied in NG demand forecasting [52–60]. In this context, ANNs were extensively used in [21,54,57,60–69] to investigate short-term NG forecasts, while in [56] different types of ANN algorithm were explored to forecast gas consumption for residential and commercial consumers in Istanbul, Turkey. ANNs were also used in [70] for the daily and weekly prediction of NG consumption of Siberia, using historical temperature and NG consumption data, in [71] for NG output prediction of USA until 2020, as well as in [72] for the prediction of NG consumption and production in China from 2008 to 2015, applying the grey theory along with NNs. Furthermore, a combination of ANNs was applied in [54,55] for the prediction of NG consumption at a citywide distribution level.

More recent studies presented numerous techniques for NG demand forecasting, including computational intelligence-based models (ANNs), fuzzy logic and support vector machines [12,73–76]. In this sense, a combination of recurrent neural network and linear regression model was used in [77] to generate forecasts for future gas demand, whereas a multilayered perceptron (MLP) neural network was deployed in [78] to estimate the next day gas consumption. A day-ahead forecast was also examined in [79] by developing a functional autoregressive model with exogenous variables (FARX). Moreover, machine learning tools such as multiple linear regression (MLR), ANN and support vector regression (SVR) were devised in [80] to project NG consumption in the province of Istanbul, as well as in [77] to forecast the residential NG demand in the city of Ljubljana, Slovenia. When considering AI methods, self-adapting intelligent grey models were also deployed for forecasting NG demand, as in [81,82].

Regarding neural network algorithms, the multilayer perceptron and the radial basis function network with different activation functions were trained and tested in [21,63,65], while the authors of [66] used a multilayer perceptron algorithm for neural network and compared this model with two time series models. In addition, Taspinar et al. explored daily gas consumption forecasting through different methods including the seasonal autoregressive integrated moving average model with exogenous inputs (SARIMAX), multi-layer perceptron ANN (ANN-MLP), ANN with radial basis functions (ANN-RBF), and multivariate ordinary least squares (OLS) [63]. Different sets of AI methods were also implemented in the following studies regarding NG consumption. ANNs with linear regression models were used in [83] for daily prediction, while ANN and Fuzzy ANN models were investigated in [84] regarding consumption in a certain region of Poland. Finally, it is worth mentioning that there has been similar research in [12], which proposed the hybrid wavelet-ANFIS/NN model to compute day-ahead forecasts for 40 distribution nodes in the national NG system of Greece.

#### *1.2. Related Work on ANFIS in Energy Consumption Forecasting*

There are plenty of works in the relevant literature regarding the application of the ANFIS model in forecasting energy consumption. As regards the electricity domain, ANFIS model was applied to forecast annual regional load in Taiwan [85] and annual demand in Turkey [86], showing in both cases that the results are good and the ANFIS model performed better than regression, neural network and fuzzy hybrid systems. ANFIS was also used in [87] for short-term electricity demand forecasting, using weekly electricity load data, as well as in [88], to estimate possible improvement of electricity consumption. Also, for electricity load forecasting, ANFIS was used in [89] to highlight its superiority to the ANN model, while it was furthermore applied in the field of transportation, forecasting the corresponding energy demand for the years 2010 to 2030, in the country of Jordan, revealing the efficiency of the examined model. Another study regarding the energy domain, where the ANFIS model was applied, is that of [90]. A long-term prediction of oil consumption was studied, further examining the interrelationship between oil consumption and economic growth in Turkey, for the years 2012 to 2030.

#### *1.3. Related Work on ANFIS in Natural Gas Consumption Forecasting*

Casting a view on the literature that refers to NG consumption forecast, the authors came across only one study that devised solely an ANFIS model. Specifically, ANFIS was used in [91] in order to estimate the daily NG demand in Iran, which actually used an extremely small dataset of historical data for both testing and training (December 2007–June 2008). Models trained on a small dataset tend to overfit, which results in high variance and very high error on a test set, producing inaccurate results. In this case, the predicting error decreases monotonically with the size of training set [92]. The rest of the studies dealt with approaches that combine ANFIS with other methods. For example, in [21], statistical time series analysis along with ANN and ANFIS methods were applied in order to predict weekly NG consumption in Turkey. Moreover, an ANFIS-fuzzy data envelopment analysis (FDEA) was developed in [93] for long-term NG consumption forecasting and analysis. In this study, 104 ANFIS were constructed and tested and six models were proposed to forecast annual NG consumption. The same approach was proposed in [94] for accurate gas consumption estimation in South America with noisy inputs. An ANFIS-stochastic frontier analysis (ANFIS-SFA) approach was formulated in [95] for long-term NG consumption prediction and analysis. Three patterns of the hybrid ARIMA–ANFIS model were tested in [2] to predict the annual energy consumption in Iran, using a set of data like population, GDP, export and import. Finally, a hybrid model of adaptive neuro fuzzy inference system and computer simulation for the prediction of NG consumption was developed in [96].

#### *1.4. Related Work on Fuzzy Cognitive Maps (FCMs) in Energy and Natural Gas Consumption Forecasting*

Moreover, other soft computing techniques, like evolutionary fuzzy cognitive maps (FCMs) have been applied for the modeling and prediction of time series problems. The dynamic modeling structure of FCMs inheriting the learning capabilities of recurrent neural networks works properly for modeling and time series prediction. Salmeron and Froelich further investigated the applicability of FCMs in univariate time series prediction by proposing an FCM simplification approach with the removal of nodes and weights [97]. Regarding the task of multivariate time series prediction, Froelich and Salmeron proposed a nonlinear predictive model based on an evolutionary algorithm for learning fuzzy grey cognitive maps [98], while Papageorgiou et al. [99] and Poczeta et al. [100] applied a new type of evolutionary FCM enhanced with the Structure Optimization Genetic Algorithm (SOGA) in energy for electricity load forecasting. Through the SOGA algorithm, an FCM model can automatically be constructed by taking into consideration any available historical data. A two-stage prediction model for multivariate time series prediction, based on the efficient capabilities of evolutionary FCMs and enhanced by structure optimization algorithms and ANNs, was introduced in [101]. In the first stage of the prediction model, SOGA-FCM was applied for selecting the most significant concepts and defining the relationships between them. Next, that model was fed into the second stage to define the initial features and weights of the training ANN. This generic prediction approach was applied in four common prediction problems, one of which dealt with electric power consumption.

In [102], Poczeta and Papageorgiou conducted a preliminary study on implementing FCMs with ANNs for NG prediction, showing for the first time the capabilities of evolutionary FCMs in this domain. Furthermore, the research team in [13] recently contacted a study for time series analysis devoted to NG demand prediction in three Greek cities, implementing an efficient ensemble forecasting approach through combining ANN, real coded genetic algorithm (RCGA)-FCM, SOGA-FCM, and hybrid FCM-ANN. In this research study, the advantageous features of intelligent methods through an ensemble to multivariate time series prediction in NG demand forecasting are explored.

#### *1.5. Research Gap and the Novelty of This Study*

Based on the reported literature survey and reviews [7,15], regarding the application of ANN-based and hybrid forecasting methods, a research gap has been identified in the field of NG day-ahead demand prediction. The observed gap mainly refers to the lack of model simplicity and flexibility, the insufficient exploration of certain modelling aspects, and inadequacy to cope with the inherent fuzziness in data handling. Most of these forecasting methods need a large dataset to be trained and a relatively large number of features to make accurate predictions. Furthermore, they are complex in their structure, time consuming and difficult to be used by non-experienced AI users. There has been hardly any research on successfully applying the ANFIS technique on the field of NG demand prediction, having performed a deep exploration process for determining the best model configuration, thus producing a highly accurate model with generalization capabilities.

Considering the aforementioned limitations, this work aims to fill the observed research gap and seeks to develop an easy to use, robust and flexible ANFIS model, which is at the same time fast, simple in structure and able to cope with fuzziness. More specifically, the proposed ANFIS architecture uses as model's inputs the most important and commonly used input variables according to the literature [7,15], such as day, month and daily average temperature, along with past NG consumption data. Moreover, a relatively large dataset was used for both testing and training the model, resulting in a systematic improvement of the model's predictive accuracy [92]. The current work pays great attention to the generalization of the proposed method and tries to properly evaluate the model's generalization capabilities in two ways: (i) by applying the model on a city level, where 10 different cities were properly examined, and (ii) by carrying out an exploration process, where 94 different models' configuration sets were examined for each one of the cities that participated in this research.

To sum up, the innovations offered in this paper are as follows, highlighting the contribution of this work to the research community:


#### *1.6. Aim of This Research Work*

The motivation of this work is to propose an ANFIS-based forecasting approach with generalization capabilities for short-term (day-ahead) city-level NG prediction in Greek areas. Also, a comparative analysis is conducted, applying ANNs, evolutionary FCMs and hybrid combinations of them on the same dataset to show the capabilities of the proposed ANFIS architecture.

The objectives of the present paper regarding NG demand forecasting, are briefly summarized in the following:

(a) To develop a robust ANFIS model to provide accurate short-term forecasts for a number of cities in Greece, using a relatively large dataset. At the same time, the authors perform model fine-tuning that can lead to high accuracy in most distribution points. The proposed model is characterized by high flexibility, easiness of use and low execution time requirements.


The outline of this paper is as follows. Section 2 describes the datasets of NG demand, collected for 10 Greek distribution points, as well as the proposed methodology of ANFIS for NG demand prediction using a well-defined set of evaluation performance metrics. Section 3 presents the results of the investigated ANFIS architectures. In the same section, a comparative analysis with other traditional neural networks and soft computing methods was performed for the same dataset. The discussion of results, which is also included in Section 3, presents the main outcomes of a meticulous ANFIS exploration analysis, along with ANFIS advantageous features. These are compared with ANNs and FCMs, and their overall contribution in NG forecasting is presented. Section 4 summarizes the paper, presenting future challenges in energy demand forecasting and highlighting further research directions.

#### **2. Materials and Methods**

This study aims to develop an ANFIS architecture capable of forecasting short-term NG consumption demand of the 10 main cities in Greece using the dataset that was provided by the Hellenic Gas Transmission System Operator S.A. (DESFA) [103]. The developed ANFIS approach deployed the aforementioned dataset along with other variables like the average daily temperature data for all the examined cities to accomplish forecasting. The results produced were further compared with those calculated by ANN and other soft computing techniques like FCM and hybrid-ANN to prove the prediction performance of the ANFIS prediction tool. Details on the dataset and its features, as well as the proposed methodology, are provided below. MATLAB M-file environment version 9.3.0.71 (R2017b) was used to program ANFIS networks and develop ANFIS models.

#### *2.1. Dataset*

The dataset covers ten different prediction datasets of historical data referring to ten cities all over Greece (Alexandroupoli, Athens, Drama, Karditsa, Larissa, Markopoulo, Serres, Thessaloniki, Trikala and Volos) and was linked to the values of gas demand for eight (8) previous years, in total. It should be mentioned that the time period for each dataset (city) was not the same in duration and did not correspond to the same years of data with all the other datasets collected. Table 1 depicts the duration in years that is linked to each dataset collected and used in this case study. The historical datasets for 15 Greek cities were initially provided by the NG grid company of Greece, DESFA, which is responsible for the operation, management, exploitation and development of the Greek NG system and its interconnections. However, the authors, after thoroughly reviewing the available datasets, decided to include only 10 out of 15 cities in their case study, since these datasets contained less outliers and missing values than the other 5 datasets that were finally rejected, for data consistency purposes. For the datasets that were finally included in this work, a preliminary preprocessing phase was performed, where the insignificant outliers were removed, and any missing values were substituted with the average real value of the previous two days demand. The real data that were used for ANFIS modeling, performance evaluation and comparison with other popular forecasting methods were then split into training and testing samples. For all cities, the last year of each dataset (from November 2017 up to October 2018) was devoted to testing, whereas the rest of the years were used for training the developed ANFIS model.


**Table 1.** Time period referred to in each time-series dataset for all cities.

In order to properly forecast day-ahead NG consumption demand of Greece, the proper number and type of input parameters should be selected. So, five factors were carefully considered as input parameters and the amount of one-day-ahead NG consumption demand of each distribution point was the output parameter. The prediction model was based on observations of past NG consumption, weather data, and calendar indicators, which are all among the most important input variables for prediction of NG consumption [15]. In particular, the dataset contains historical data of NG consumption of each city's distribution point, the daily average temperature of the area in Celsius degrees, a month indicator and a day indicator. As regards the previous NG consumption data, these are linked to two different input variables: demand of a day before and current day demand. The temperature data are obtained by the nearest to the distribution gas point meteorological station. Concerning the calendar indicators (month and day), they need to undergo certain data form preprocessing before their use. Specifically, two different input indicators need to be considered for each one of the two variables. We define k = 1,2, ... , 12 as the month index (1 January, 2 February, ... , 12 December) and l = 1,2, ... , 7 as the day index (1 Monday, 2 Tuesday, ... 7 Sunday). Following the coding procedure as presented in [104], the index for the month is scaled to the range [1/12, 1] in which the months of the year from January to December take successive values of the scaled index. That is, January has the value of 1/12 and December the value of 1. Similarly, the days of the week take successive values in the scaled range [1/7, 1], in which Monday and Sunday take the values of 1/7 and 1, respectively. All these parameters constituting the actual recorded data are briefly presented in Table 2.


**Table 2.** Input and output parameters.

All data that compose the investigated dataset underwent a normalization process. This was necessary because all entries needed to have the same limited range of values so the model produces meaningful results [105].

The algorithm that was used for data normalization is the Min-Max, which scales the values of the dataset linearly over a specific range. As described in previous works [13,105], each variable was normalized to [0,1] before the forecasting model was applied. The normalized variable took its original value when the testing phase was implemented. Data normalization was carried out mathematically, as follows:

$$\mathbf{x}\_i^{(new)} = \frac{\mathbf{x}\_i - \mathbf{x}^{(min)}}{\mathbf{x}^{(max)} - \mathbf{x}^{(min)}}, \forall i = 1, 2, \dots, N \tag{1}$$

where *x*(*new*) is the normalized value of the variable *x*, and *x*(*min*) and *x*(*max*) are, respectively, the minimum and maximum values of the concerned variable *x*.

#### *2.2. Methods*

#### 2.2.1. Adaptive Neuro-Fuzzy Inference System (ANFIS)

The adaptive neuro-fuzzy inference system (ANFIS) uses an architecture that is based on both ANN and fuzzy logic principles and takes advantage of the benefits of both in a single framework. It can be described by the fuzzy "IF-THEN" rules from the Takagi and Sugeno (TS) type [106] as follows:

$$\begin{aligned} R\_i: \text{if } \mathbf{x}\_1 = A\_{i,1} \text{ and } \dots \text{and } \mathbf{x}\_k = A\_{i,k} \\ \text{then } y\_i = b\_{i,0} + b\_{i,1}\mathbf{x}\_1 + b\_{i,2}\mathbf{x}\_2 + \dots + b\_{i,k}\mathbf{x}\_k \end{aligned} \tag{2}$$

where *Ai*,*<sup>k</sup>* is the membership function associated with input variables *xk* and *n* is the number of inputs.

A typical ANFIS network is a five-layer structure consisting of the fuzzy layer, the product layer, the normalized layer, the de-fuzzy layer and the total output layer [3,107,108], as depicted in Figure 1.

**Figure 1.** The TSK ANFIS architecture [107].

In the first layer, every node *i* represents a linguistic label and is described by the following membership function, as given in Equation (3).

$$A\_{i,k}(\mathbf{x}\_r) = e^{-\left(\frac{\mathbf{z}\_r - v\_{i,k}}{\sigma\_{l,k}}\right)^2}, \text{ for } r = 1, 2, \dots, i \tag{3}$$

where *Ai*,*<sup>k</sup>* is the membership function which is considered to be Gaussian and is described by the center ν and the spread σ.

In the second layer, the firing strength of the rule is computed using multiplicative operator, as presented in Equation (4). Firing strength is the weight degree of the IF-THEN rule and determines the shape of the output function for that rule.

$$w\_i = \prod\_{k=1}^{n} A\_{i,k}(x\_k) \tag{4}$$

In the third layer, the *i*-th node calculates the ratio of the *i*-th rule's firing strength to the sum of the firing strength of all rules. This is the normalization layer which normalizes the strength of all rules and the output of each node is given by Equation (5).

$$
\overline{w}\_{i} = \frac{w\_{i}}{\sum\_{i=1} w\_{i}} \tag{5}
$$

In the fourth layer, each node is an adaptive node with a function given by Equation (6). In this layer, each node calculates a linear function where its coefficients are adapted by using the error function of the multilayer feed-forward neural network.

$$
\overline{y}\_i = \overline{w}\_i y\_i \tag{6}
$$

In the fifth layer, there is only a fixed node indicated as the sum of the net outputs of the nodes in Layer 4. It computes the overall output as the sum of all incoming inputs and is expressed by Equation (7).

$$y = \sum\_{i} \overline{y}\_i \tag{7}$$

ANFIS uses a hybrid learning algorithm to train the model. The back-propagation algorithm is used to train the parameters in Layer 1, whereas a variation of least-squares approximation or back-propagation algorithm is used for training the parameters of the fourth layer [108,109].

#### 2.2.2. Proposed ANFIS Architecture Applied in Natural Gas Consumption Forecasting

In order to develop an efficient ANFIS model for NG demand forecasting, the authors needed to follow a certain process regarding the design of model's architecture as well as an exploration process that will properly configure the input and training parameters of the examined model. Priority was given to the definition of the FIS architecture before the training of the network [110]. Among various fuzzy inference system (FIS) models, the Sugeno fuzzy model is the most widely used because of its higher interpretability and computational ability, that includes embedded optimal and adaptive techniques [111]. In order to create a fuzzy rule, the input space needs first to be divided. Two methods are used to divide space, comprised by input variables: the grid partitioning method and the subtractive clustering method. The main difference between these two functions refers to the way the partition of the input space is created.

In grid partitioning [109], the input space is divided into a grid-like structure without overlapping parts. Grid partitioning performs partitioning of the input space using all possible combinations of membership functions of each variable. This method is used when the number of input variables is small. For example, for 10 input variables and two membership functions for each input variable, then the input space is divided into 2<sup>10</sup> = 1024 specific areas, representing one rule for each specific area, and the total number of rules is 1024, which is a very complicated structure. Therefore, the grid partitioning method is mainly used when the number of input variables is small.

On the other hand, the subtractive clustering method divides the input space into appropriate clusters, even if the user does not specify their number. If the size of the cluster becomes small, then the number of clusters increases, thus increasing the number of fuzzy rules. A rule is created for each cluster, whereas different values for parameters, like range of influence, squash factor, accept ratio and reject ratio, need to be explored for determining an efficient architecture, which will keep the balance between the total number of ANFIS parameters and the total number of rules.

Considering the above specifications, the authors used the Grid partition option to define the FIS architecture due to its simplicity, less time-consuming performance as well as it can easily explore the number and type of membership function (MF). In this stage, the number and type of membership functions of each input variable, along with the rules and values of parameters that belong to these functions, were determined using the option of Grid partition.

When implementing an ANFIS architecture, researchers should have in mind that there is one main restriction: the number of input variables. When these are more than five, then the number of the IF-THEN rules and the computational time also increase, hindering ANFIS to model output with respect to inputs [110]. Thus, in this study, five variables were chosen as input parameters, i.e., month, day, temperature, demand of a day before and demand of current day. As described above, a day-ahead consumption demand was selected as the output variable whose value can be produced by choosing between the option of linear or constant type of MF.

Finding the most efficient ANFIS architecture is a demanding task and entails a rigorous exploration process. Since our concern focuses on the increment of network's accuracy and decrement of the errors, five necessary configurations should be considered in this direction: (i) the number of membership functions, (ii) types of MF (triangular, trapezoidal, bell-shaped, Gaussian and sigmoid), (iii) types of output MF (constant or linear), (iv) optimization methods (hybrid or back propagation) and (v) the

number of epochs [112]. For the convenience of readers, these steps are visually represented in the flowchart in Figure 2.

**Figure 2.** Flowchart of the proposed adaptive neuro-fuzzy inference system (ANFIS) methodology.

The aforementioned set of configurations needs to be deployed in order to generate FIS and next to train the ANFIS model. Accordingly, the dataset that included the five input variables (i.e., month, day, temperature, demand of a day before, demand of current day) was selected to determine the only output (day-ahead demand). Initially, the training dataset was loaded in the ANFIS tool, as shown in Figure 3a. The next step was the design of the neuro-fuzzy model using the option "Generate FIS". The Grid partition option was selected according to the description above (see Figure 3a). These two settings, concerning the fuzzy input variables along with their membership functions, are the most important parts to design the ANFIS. An example of selecting the number and type of MFs is illustrated in Figure 3b.

**Figure 3.** Screenshots regarding (**a**) the training dataset configuration, (**b**) the number and type of MF configuration.

The number and type of membership function were assigned to the input parameters following the trial-and-error approach. The different types of MF that are offered by the MATLAB ANFIS editor include the triangular, trapezoidal, generalized bell (Gbell), Gaussian curve, Gaussian combination, difference between two sigmoid functions and product of two sigmoid functions (see Figure 3b). Regarding the type of output MFs, in the Sugeno-type fuzzy system, there are two options: a constant-type conclusion or a linear-type conclusion function. In the case of linear function, the output *y* is defined as:

$$y = k\_0 + k\_1 \ast x\_1 + k\_2 \ast x\_2 + \dots + k\_n \ast x\_n \tag{8}$$

where *x*1, *x*2, ... , *xn* are the *n* inputs. In this case, ANFIS needs to define *k*0, *k*1, *k*<sup>2</sup> up to *kn*, and it is very time consuming to efficiently calculate the outputs when a large number of parameters are considered. On the other hand, when a constant MF is selected, the algorithm needs to define only one parameter to provide a reliable forecasted value. Thus, the computational time is really low.

The selected configuration also includes the hybrid optimization method, while the number of epochs selected to train the model was between 10 and 50. The hybrid optimization method uses the back propagation learning algorithm for parameters associated with input MF and the least-square estimation algorithm for parameters associated with output MF; thus, it was selected as the most proper one [113]. Various sets of ANFIS configurations are presented in Table 3, regarding different sets of number for MFs, as considered by the authors of this work.

**Table 3.** Different configurations of the selected ANFIS architectures regarding constant outputmembership function (MF).


Regarding the output MFs, constant and linear MF were accordingly investigated after certain numbers of experiments conducted. From these experiments and for the linear output, it was observed that the number of rules increases significantly, as well as the computational time, even in the case of problems with a small number of inputs (see Table A1 in Appendix A). Thus, the linear type was not considered as an appropriate parameter of output MF since it is extremely time consuming. In this context, the trial-and-error approach was followed for the selection of the input-output type of MFs. Figure 4 illustrates an indicative ANFIS model, which was constructed with the following configuration set: 3-3-3-2-2, gbell MF, constant output MF, 10 epochs, hybrid.

**Figure 4.** Screenshot of the constructed ANFIS model.

#### 2.2.3. Testing and Evaluation

The testing process for the ANFIS model was accomplished by using the testing data, which were completely unknown to the model. The predictor makes predictions on each day and finally compares the calculated predicted value with the real value. For example, considering the city of Volos, the predicted values that are illustrated in red in Figure 5 are compared with the real values (in blue color).


**Figure 5.** Screenshot of the testing data configuration.

In order to evaluate the prediction of NG demand, five well known and commonly used statistical indicators were introduced, i.e., mean square error (MSE), root mean square error (RMSE), mean absolute error (MAE), mean absolute percentage error (MAPE) and coefficient of determination (R2). The mathematical equations of the statistical indicators are described below.

1. Mean squared error:

$$\text{MSE} = \frac{1}{T} \sum\_{t=1}^{T} (Z(t) - X(t))^2 \tag{9}$$

*Energies* **2020**, *13*, 2317

2. Root mean squared error:

$$\text{RMSE} = \sqrt{\text{MSE}} \tag{10}$$

3. Mean absolute error:

$$\text{MAE} = \frac{1}{T} \sum\_{t=1}^{T} \left| Z(t) - X(t) \right| \tag{11}$$

4. Mean absolute percentage error:

$$\text{MAPE} = \frac{1}{T} \sum\_{t=1}^{T} \left| \frac{Z(t) - X(t)}{Z(t)} \right| \tag{12}$$

#### 5. Coefficient of determination:

$$\mathcal{R} = \frac{T\sum\_{t=1}^{T} Z(t) \cdot X(t) - \left(\sum\_{t=1}^{T} Z(t)\right) \left(\sum\_{t=1}^{T} X(t)\right)}{\sqrt{T\sum\_{t=1}^{T} \left(Z(t)\right)^2 - \left(\sum\_{t=1}^{T} Z(t)\right)^2} \cdot \sqrt{T\sum\_{t=1}^{T} \left(X(t)\right)^2 - \left(\sum\_{t=1}^{T} X(t)\right)^2}} \tag{13}$$

where *X*(*t*) is the forecasted value of the NG at the *t*-th iteration, and *Z*(*t*) is the actual value of the NG at the *t*-th iteration, *t* = 1, ... , *T*, where *T* is the number of testing records.

Higher values of R2, i.e., closer to 1, mean better model performance and the regression line fits the data well. A coefficient of determination value of 1.0 points out that the regression curve fits the data perfectly.

#### **3. Results**

This section presents the exploration analysis results for the various ANFIS architectures as proposed in Section 2.2.2. Considering the steps proposed in Section 2.2, the initial dataset is split into training and testing. During the training process, the ANFIS model is designed for each one of the suggested configurations. After the training process of the ANFIS finishes, NG consumption demands for the next day (one day ahead prediction) are calculated from the generated FIS. The NG consumption results only for the city of Athens (as an indicative example) regarding all configurations tested, are presented in Table 4, whereas Table 5 gathers the best three results of NG consumption demand obtained from ANFIS for each of the 10 cities. To evaluate the performance of the models, the statistical indicator mean absolute percent error (MAPE) was used [91]. MAPE is a relative measurement, independent of scale, and it is the most common performance metric in time series forecasting, due to being reliable and valid [21].

**Table 4.** Testing results for the examined ANFIS architectures concerning the city of Athens.


**Table 4.** *Cont.*


**Table 5.** Testing results for the best three ANFIS architectures of each city based on MAPE value.


The corresponding graphical representation of the results regarding the best three out of 47 total ANFIS architectures for the city of Athens is illustrated in Figure 6. Also, the best ANFIS model for each city can be found in Table 6, which provides the most reliable ANFIS architecture results for each city. All the results have been previously ranked, based on the minimum value of MAPE and subsequently the minimum values of MSE, RMSE and MAE. The priority was given to MAPE as one of the most crucial evaluation metrics, according to the literature [91,96], which was used in this study to compare various models obtained from ANFIS and other soft computing and neural networks methods. As a relative and easy to interpret measurement, MAPE is reliable, valid and independent of scale. The smaller the values of MAPE are, the closer the forecasted values are to the actual values.

**Figure 6.** Forecasting results for the best three ANFIS architectures for Athens.


**Table 6.** The best ANFIS models for all the cities under investigation (10 epochs and hybrid optimization).

As illustrated in Table 6, ANFIS models appear to perform best mostly when triangular MFs are used for the input variables: three MFs for the first three input variables (month, day of week and mean temperature) and two or three MFs for the other two input variables (daily demand for current day and one day before). Also, constant MFs are selected for the output variable and hybrid optimization method. The graphical representation of the best ANFIS models for each city is illustrated in Figure 7.

**Figure 7.** *Cont.*

**Figure 7.** Forecasting results for three cities considering the best ANFIS method. (**a**) Testing Alexandroupoli, (**b**) testing Athens, (**c**) testing Drama, (**d**) testing Karditsa, (**e**) testing Larissa, (**f**) testing Markopoulo, (**g**) testing Serres, (**h**) testing Thessaloniki, (**i**) testing Trikala, (**j**) testing Volos.

It is worth mentioning that all three most efficient ANFIS architectures with respect to MAPE values have triangular or gaussian MFs and 2-2-2-2-2 or 3-3-3-2-2 number of input MFs, whereas the output MF is constant and the learning algorithm is hybrid. In addition, the application of other MFs combinations does not seem to give results that could be on top of the list. Due to the limitation of the total number of parameters that should not exceed the number of training data pairs, the number of MFs was chosen based on the number of input parameters. Figure 8 shows the exponential increase in the number of rules when the number of MFs increases, whereas Table A2 in Appendix A gathers the time and number of rules for all the proposed ANFIS configurations.

**Figure 8.** Number of rules for various ANFIS architectures.

#### *3.1. Comparison with ANNs, FCMs and Hybrid FCM-ANN*

To further investigate the performance of the proposed ANFIS architectures, an extensive comparative analysis between the state of the art ANNs, soft computing methods of FCMs and their hybrid combination of FCMs with ANNs was performed.

The architecture of the analyzed ANN was a multilayer feed forward network with an input layer containing five inputs (month, day, temperature, demand of a day before, current demand), a hidden layer with 10 neurons, and an output layer with one output (a day-ahead demand prediction). The authors used the sigmoidal activation function in all layers and implemented Levenberg–Marquardt algorithm to train the network.

The soft computing method of fuzzy cognitive map with evolutionary learning capabilities, such as the real-coded genetic algorithm (RCGA-FCM) and structure optimization genetic algorithm (SOGA-FCM) [100], were used for time series modeling and prediction of day-ahead NG energy demand. For FCM learning, we implemented RCGA-FCM and SOGA-FCM. A short description on the applied evolutionary-based FCM approaches is given in Appendix B. The implementations of FCMs differ from ANFIS, even though they both belong to the soft computing family.

In this research study, we used the dynamic model type (Equation (B1)) which is found in Appendix B, with sigmoidal transformation function. FCMs learned with the use of RCGA and SOGA algorithm contain five concepts (month, day, temperature, demand of a day before, current demand) [114].

The applied hybrid approach for time series prediction is based on FCMs and ANNs and was previously proposed in [18,19]. It allows us to select the most significant concepts for FCM using SOGA. These concepts are used as the inputs for ANN. In the hybrid approach, we used artificial neural networks with an input layer with five inputs selected by the SOGA-FCM approach, a hidden layer with 10 neurons and an output layer with one output (one day-ahead demand prediction). Sigmoidal activation function and Levenberg–Marquardt learning algorithm were used. All the simulations for FCMs and hybrid FCM-ANN configurations were performed with the software tool ISEMK [115] which has been developed for time series forecasting purposes. An analytical description of FCM-based models and hybrid FCM-ANN can be found in [13,101,102,116], whereas they are used in this work only for comparison purposes.

In what follows, Table 7 gathers the results of the explored ANN and soft computing models, which are straightforward compared with our best performed ANFIS configuration, for each one out of the 10 cities, suggested in this research work. In Figure 9, three indicative graphs of the cities Alexandroupoli, Athens, and Drama are illustrated regarding the predicted values of NG demand for

all the best proposed architectures. In Appendix A, the corresponding graphs of the rest of the cities are presented (see Figure A1).


**Table 7.** Comparison results among the artificial neural network (ANN), fuzzy cognitive map (FCM), hybrid FCM-ANN and best ANFIS architectures of each city.

**Figure 9.** Comparison of forecasting results for each city considering all examined methods. (**a**) Testing for Alexandroupoli, (**b**) Testing for Athens, (**c**) Testing for Drama.

For a deeper analysis of the examined architectures (ANFIS, ANN, FCM, hybrid FCM-ANN), the authors report on further details regarding the parameters of each model used in this study. The ANN and FCM models were previously applied for NG demand prediction in several research works, such as those in [13,101,102,116]. The models were sufficiently described and the hyperparameters were properly configured to offer optimum performance of the investigated FCM models.

Table 8 depicts the optimum parameters for all cities considering the neural and FCM evolutionary methods (ANN, RCGA-FCM, SOGA-FCM, Hybrid), compared with the proposed best performed ANFIS. The average running time is also presented in Table 8, which was calculated for each soft computing architecture for all models. It is worth mentioning that we have conducted a rigorous exploratory analysis for all the investigated neuro-fuzzy, soft computing techniques and ANNs, with different parameters, for training and model optimization, to reach the highest prediction accuracy with respect to the evaluation metrics.


**Table 8.** Parameters and average running time for each architecture.

#### *3.2. Discussion of Results*

In this work, several ANFIS architectures were investigated, with respect to all the variables that were carefully determined in the developed model as reported in Section 2.1 and after different sets of model configurations were tested. However, only one ANFIS architecture reached the optimum performance, in terms of forecasting accuracy, considering the minimum value of MAPE and subsequently the minimum values of MSE, RMSE and MAE values produced. In particular, it emerged that the optimum ANFIS configuration is 2-2-2-2-2 with triangular MFs for input variables, which produces the most simple (concerning the number of rules), fast (see Table A2 in Appendix A) and accurate model for this energy forecasting problem. In general, it is observed that the best results are produced from the combination of triangular or gaussian MFs regarding the input variables, and the constant MFs regarding the output layers.

To further discuss the results produced and to show the effectiveness of the proposed forecasting methodology of ANFIS, the authors conducted a comparative analysis regarding the forecasting performance between the proposed technique, and other ANN and soft computing methods too, such as FCM, which were reported in the literature and have already been applied in the specific domain. The MAPE criterion was used to compare various models from ANN, evolutionary FCM, Hybrid FCM-ANN and ANFIS. The results are given in Table 7. For example, the MAPE values of RCGA-FCM, SOGA-FCM, ANN, hybrid and ANFIS models for the city of Alexandroupoli were calculated as 17.62%, 17.17%, 16.11%, 14.30% and 10.11%, respectively. The smaller the values of MAPE are, the closer to the actual values the forecasted values are. The best result was obtained from the ANFIS model. The respective figures in the text and in Appendix A have been updated with the new prediction values.

Considering the same dataset linked to only three cities (Athens, Thessaloniki, and Larissa) out of the ten that participated in our study, a day-ahead NG consumption prediction was investigated in [117], applying ANN and LSTM approaches and in [102], implementing the SOGA-FCM method and a hybrid combination of it. Furthermore, an ensemble FCM prediction methodology concerning the same dataset was presented in [13], in which a recent soft computing technique for time series forecasting, using evolutionary fuzzy cognitive maps and their ensemble combination was compared to ANNs, as benchmark forecasting methods. These methods and their results in terms of MSE and MAE values for three benchmark cities are all gathered in the following table and certain results can be concluded. The main reason for selecting the statistical indicators MSE and MAE in the following figure is to accomplish a straightforward comparison with the results published in previous works.

In Figure 10, it can be noted that all methods achieve high accuracy in NG consumption predictions, using the same dataset. The best ANFIS approach seems to excel over the ensemble and hybrid methods. Consequently, the proposed ANFIS architecture, which handles the fuzziness of data more efficiently, outperforms all the other examined methods in most cases, with a rather remarkable difference. ANFIS is less time consuming and more flexible than ANN, and as it employs fuzzy rules and membership functions incorporating with real-world systems, it can be used as alternate method to ANN forecasting.

**Figure 10.** Comparison of results between machine learning and soft computing methods for three benchmark cities.

It is observed that the proposed method exhibits better or similar performance to other well-known ANN, FCM or hybrid FCM-ANN architectures for the ten cities under investigation. The produced results highlight the significance and superiority of neuro-fuzzy methods over the other examined methods in terms of prediction accuracy, when they deal with time series forecasting problems in energy. This is in accordance with the main advantageous features of ANFIS models, which are their ability to capture the nonlinear structure of a process, their adaptation capability, and fast training characteristics. As reported in the literature, ANFIS models are able to cope with the uncertainty and fuzziness that characterize the energy domain [118,119] when other intelligent methods cannot tackle them.

The main outcomes of this study can be summarized as follows:


#### **4. Conclusions**

This study proposes the ANFIS method to predict short-term demand of NG consumption. This approach is applied on the Greek territory and uses 10 different datasets provided by DESFA, that regard previous energy consumption historical data, for ten main cities. To decide the model's proper architecture, the authors follow an exploration process regarding the best configuration of input and training parameters. The best ANFIS model is then compared to other well-known ANN and soft computing models that are commonly used for energy demand prediction purposes. The ANFIS method demonstrates significant performance in the field of energy demand prediction, outweighing the traditional ANN and FCM architectures. In addition, the running time of the proposed architecture is much less than those of other examined models, making it the right decision for day-ahead demand forecasting of NG. The findings of this study reveal that the highest forecasting accuracy emerged when the same model configuration was used for most of the cities, highlighting the generalization capabilities of the proposed architecture.

This work can be widely used in short-term demand forecasting for other countries too, with the same or similar input parameters, and can be also useful especially for distribution operators, providing them with the ability to make long-term planning decisions and apply the correct strategic policies in this direction. Following the literature, the ANFIS approach can be applied in various other domains such as medicine, environmental modelling, various energy systems, like solar and wind, as well as other engineering applications. As can be seen, the ANFIS application area is wide, and as regards the energy sector, this method finds great applicability due to its high prediction accuracy, robustness, and easiness to use.

The results show that the proposed algorithm, which was proven to be efficient, fast and robust, can be adopted by regulatory authorities and decision makers to perform rigorous forecasting of natural gas demand for the respective case cities and other cities in Greece too. The investigated approach is an accurate estimation method as it makes efficient short-term predictions in natural gas demand, showing minor deviations between the real and the predicting values. Since short-term natural gas forecasting is mostly used for the timely reservation of transport, storage capacity optimization, timely purchase of natural gas deliveries and capacity allocation, this method becomes critical to determine the energy policy for Greece and the wider area too, having overall a positive impact in natural gas consumption.

Future work is oriented in developing more advanced neuro-fuzzy models providing explainability and transparency in prediction tasks in diverse research domains, in order to evaluate the generalization capabilities of this approach. Furthermore, new forecast combination architectures of efficient deep learning and regularized recurrent neural networks for time series modelling and prediction in the energy sector will be investigated.

**Author Contributions:** Conceptualization, E.I.P. and E.I.P.; methodology, E.I.P.; software, E.I.P., K.P. (Konstantinos Papageorgiou) and K.P. (Katarzyna Poczeta); validation, E.I.P., E.I.P., K.P. (Konstantinos Papageorgiou), K.P. (Katarzyna Poczeta), D.B. and G.S.; formal analysis, E.I.P. and E.I.P.; investigation, E.I.P.; resources, E.I.P.; data curation, E.I.P. and K.P. (Konstantinos Papageorgiou), K.P. (Katarzyna Poczeta); writing—original draft preparation, E.I.P.; writing—review and editing, K.P. (Katarzyna Poczeta), E.I.P., D.B. and G.S.; visualization, E.I.P.; supervision, E.I.P. and G.S.; project administration, E.I.P. and D.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

#### **Appendix A**


**Table A1.** Different configurations of the selected ANFIS architectures regarding linear output MF.



**Figure A1.** *Cont.*

**Figure A1.** 76 *Cont.*

**Figure A1.** Comparison of forecasting results for each city considering all examined methods. (**a**) Testing for Karditsa, (**b**) testing for Larissa, (**c**) testing for Markopoulo, (**d**) testing for Serres, (**e**) testing for Thessaloniki, (**f**) testing for Trikala, (**g**) testing for Volos.

#### **Appendix B**

#### *Appendix B.1. Fuzzy Cognitive Maps*

Fuzzy cognitive maps (FCMs) are an effective tool for modeling and predicting time series. The structure of the FCM model is based on a directed graph, the nodes of which denote concepts significant for the analyzed problem, and the links are the causal relationships. Values of concepts can change over time according to the adopted dynamics model, for example the nonlinear dynamics model:

$$\mathbf{X}\_{i}(t+1) = F\left\{\mathbf{X}\_{i}(t) + \sum\_{\substack{j=1\\j\neq i}}^{n} \mathbf{X}\_{j}(t) \cdot w\_{j,i} \right\} \tag{A1}$$

where *Xi*(*t*) is the value of the *i*-th concept at the *t*-th iteration, *wj,i* is the weight of the causal relationship between concepts *Xj* and *Xi* taking values from the range [−1,1], *t* is discrete time, *i,j* = 1, 2, ... , *n, n* is the number of concepts, and *F* is the transformation function normalizing the factor values to the range [0,1] or [−1,1]. Fuzzy cognitive maps can be constructed based on expert knowledge or with the use of machine learning algorithms. The aim of fuzzy cognitive map learning is to determine the weights of the causal relationships between concepts on the basis of available time series.

An effective method for fuzzy cognitive map learning is the real-coded genetic algorithm (RCGA) [100]. RCGA defines each individual in the population based on a floating-point vector containing the causal relationships. Each individual is decoded into a candidate map and evaluated with the use of proper fitness function. We used the following fitness function:

$$fitness\_p(error\_l) = \frac{1}{a \cdot error\_l + 1} \tag{A2}$$

where *a* is a parameter, l is the number of generations, *l* = 1, ... *, L*, *L* is the maximum number of generations, *p* is the number of individuals, *p* = 1, ... , *P*, *P* is the population size, and *erorrl* is the learning error that can be in the following form:

$$error\_l = \frac{1}{T} \sum\_{t=1}^{T} \left( Z(t) - X(t) \right)^2 \tag{A3}$$

where *X*(*t*) is the predicted value of the decision concept at the *t*-th iteration, *Z*(*t*) is the real normalized value of the decision concept at the *t*-th iteration, *t* = 1, ... , *T*, and *T* is the number of learning records.

Another way of learning fuzzy cognitive maps is the structure optimization genetic algorithm (SOGA) [100,120]. This allows one to simplify the structure of the FCM model by selecting the most significant concepts and causal relationships during the learning process. In this approach, the fitness function is based on the modified learning error including an additional penalty for highly complexity of the candidate fuzzy cognitive map, understood as a large number of concepts and non-zero relationships, described as follows:

$$error\_{\parallel} = error\_{\parallel} + b\_1 \frac{n\_{\parallel}}{n^2} error\_{\parallel} + b\_2 \frac{n\_{\odot}}{n} error\_{\parallel} \tag{A4}$$

where *b*1, *b*<sup>2</sup> are the learning parameters, *nc* is the number of the concepts in the candidate FCM model, *nr* is the number of the non-zero relationships between concepts, *n* is the number of all possible concepts, and *erorrl* is the learning error type (Equation (A3)).

In this paper, we also used the hybrid approach for time series prediction based on fuzzy cognitive maps with the structure optimization genetic algorithm and artificial neural networks [101]. In the first stage of this approach, the most important concepts are selected with the use of FCMs and the SOGA algorithm. In the second stage, these concepts are used as inputs for the artificial neural network in order to increase the prediction accuracy. The above algorithms have been implemented in the developed ISEMK system [101].

#### **References**


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

## *Article* **Dry Above Ground Biomass for a Soybean Crop Using an Empirical Model in Greece**

#### **Christos Vamvakoulas, Stavros Alexandris and Ioannis Argyrokastritis \***

Laboratory of Agricultural Hydraulics, Department of Natural Resources Development and Agricultural Engineering, Agricultural University of Athens, 75 Iera Odos str., 11855 Athens, Greece; chrisvamva@aua.gr (C.V.); stalex@aua.gr (S.A.)

**\*** Correspondence: jarg@aua.gr

Received: 1 November 2019; Accepted: 24 December 2019; Published: 1 January 2020

**Abstract:** A new empirical equation for the estimation of daily dry above ground biomass (D-AGB) for a hybrid of soybean (*Glycine max* L.) is proposed. This equation requires data for three crop dependent parameters; leaf area index, plant height, and cumulative crop evapotranspiration. Bilinear surface regression analysis was used in order to estimate the factors entering in the empirical model. For the calibration of the proposed model, data yielded from a well-watered soybean crop for the year 2015, in the experimental field (0.1 ha) of the agricultural University of Athens, were used as a reference. Verification of the validity of the model was obtained by using data from a 2014 cultivation period for well-watered soybean cultivation (100% of crop evapotranspiration water treatment), as well as data from three irrigation treatments (75%, 50%, 25% of crop evapotranspiration) for two cultivation periods (2014–2015). The proposed method for the estimation of D-AGB may be proven as a useful tool for estimations without using destructive sampling.

**Keywords:** dry above ground biomass; soybean; empirical models; bilinear regression analysis

#### **1. Introduction**

Nowadays, agronomists and irrigation experts use crop productivity models for the simulation and prediction of dry above ground biomass (D-AGB). Complex crop growth models (CGM) require a large number of input parameters, usually not available from ideal sites, which leads to significant and systematic cumulative errors in determining crop yield and above ground biomass.

In this study we tried to present a simple model using geostatistical methods in a simple form. A similar approach is the responsive surface methodology (RSM), which has been used by researchers as an optimization technique with a wide range of applications, such as in dehydration operations of selected agricultural crops [1] and for biodiesel production [2–6]. Also, in the optimization of tomato yield [7], to increase resistant starch in vermicelli [8] and optimization of the plough working surface [9]. Other fields of use for this optimization method are in the use of agricultural waste [10] and for harvesting losses in corn seed [11].

Generally, empirical equations are a tool for local estimations of attributes without many parameters as inputs. In the past, algorithms for the creation of such equations have been used for the estimation of reference evapotranspiration (ET0) [12,13] and for the estimation of crop evapotranspiration [14].

Also, empirical bilinear regression equations have been used for the prediction of human and rat tissue [15], while multiple regression analysis has been applied for un-mixing of surface temperature data in an urban environment [16].

In this paper, a model for the daily estimation of D-AGB for a soybean hybrid (PR91M10) in central Greece is formulated. The model has been parameterized by experimental observations on the soybean crop. Also, the model was examined for water stressed and non-water stressed plants, under field conditions. The final equation obtained is based on leaf area index (LAI), plant height (*h*c), and cumulative crop evapotranspiration (cumETc).

#### **2. Methodology**

The experiment was performed in the experimental field of the agricultural University of Athens in Aliartos plain (38◦23 40 N, 23◦05 08 E and 95 m altitude), during 2014 and 2015 cultivation periods. Data from an experimental design with irrigation treatments (100%, 75%, 50%, 25% of crop evapotranspiration, respectively) as the main plot in a randomized complete block design, with four replications, were used. The plot size of each irrigation treatment was 3 m × 12 m and the spacing between each main plot was 3 m in order to minimize water movement among treatments. The experimental plots were 3 m × 6 m and consisted of 5 rows 0.75 m apart. PR91M10 is a highly productive variety of the early maturity group (00). Seeds were hand-planted, using a seeding depth of about 3 cm, on 30 May 2014 (Julian day, JD: 150) and on 31 May 2015 (Julian day, JD: 151), respectively. Treatment plots consisted of 5 rows planted, 75 cm apart, with 4–5 cm row spacing, and the sowing density was 33 seeds m<sup>−</sup>2. At sowing, basal fertilization was performed in all plots and corresponded to 100 and 50 kg ha−<sup>1</sup> of P2O5 and K2O, respectively, based on soil fertility analysis. Irrigation scheduling was based on the daily water balance calculation and on results obtained by using the computer model ISAREG [17], which utilized data collected during consecutive cultivation periods from 2011 to 2015. The following input data for the ISAREG model were used:


Rainfall during the 2014 cultivation period was 46.1 mm and 176.7 mm for the 2015 cultivation period. The ground water table was at 1.2 m depth for both cultivation years. Irrigation was applied to provide 100%, 75%, 50%, and 25% of the crop evapotranspiration needs.

A surface drip irrigation system was used for irrigation. A 16 mm diameter polyethylene pipe with inline pressure compensating drippers at 33 cm intervals was placed on one side of each soybean row. The average discharge of emitters was 4.4 L/h at 0.1 MPa.

Periodically, every 7 days approximately, plant height (*h*c, cm) was measured and destructive sampling was performed by collecting three plants from the three interior rows of each plot, for leaf area index (LAI) and dry above ground biomass (D-AGB, ton/ha). Sampling was performed at 25, 34, 41, 48, 54, 60, 66, and 75 days after planting (DAP) for 2014 cultivation period and at 24, 33, 40, 47, 53, 59, 65, and 75 (DAP) for the 2015 cultivation period, respectively.

The parameterization of the model was done for the 2015 cultivation year data, because precipitation was higher than that of 2014, giving better environmental conditions for the non-water stressed plants (100% treatment). The model represents the simulation curve of the D-AGB for the first 75 days of the growing period. The last 20 days of the maturity stage are not included in the simulation curve. Surface regression analysis was used to establish the new model to simulate daily (D-AGB). The empirical model was derived by surface polynomial regression using three crop dependent parameters, measured values of leaf area index (LAI), plant height (*h*c), and cumulative crop evapotranspiration (cumETc), in a general form D-AGB = *f*(LAI,hc, cumETc). It utilizes four unknown parameters (*k*0, *k*1, *k*2, *k*3) which are determined in a three stage approach. Experimental lines for the D-AGB obtained from the destructive sampling performed are used as standard values. Calculated D-AGB values are then regressed against mean daily values of pairs of LAI and *h*<sup>c</sup> (first stage) and LAI and cumETc (second stage) in a bilinear equation of the form:

$$z = f(\mathbf{x}, \mathbf{y}) = k\_0 + k\_1 \cdot \mathbf{y} + k\_2 \cdot \mathbf{x} + k\_3 \cdot \mathbf{x} \cdot \mathbf{y}$$

*x*, *y* denoting daily values of either LAI and *h*c (cm), in the first stage of investigation, or LAI and cumETc (mm/time), in the second stage, *z* standing for D-AGB (ton/ha). As expected, the first and second stages end up with the estimation of two sets of four parameters *ai*, *bi* (*i* = 1, ... 4), as shown in the Equations (1) and (2) below:

$$C\_1 = a\_1 + a\_2 \cdot h\_\circledast + a\_3 \cdot \text{LAI} + a\_4 \cdot \text{LAI} \cdot h\_\circledast \tag{1}$$

where *a*<sup>1</sup> = −0.143, *a*<sup>2</sup> = 0.095, *a*<sup>3</sup> = −6.33, *a*<sup>4</sup> = 0.058.

$$C\_2 = b\_1 + b\_2 \cdot \text{cumET}\_c + b\_3 \cdot \text{LAI} + b\_4 \cdot \text{LAI} \cdot \text{cumET}\_c \tag{2}$$

where *b*<sup>1</sup> = −0.115, *b*<sup>2</sup> = 0.0066, *b*<sup>3</sup> = −2.4, *b*<sup>4</sup> = 0.0129.

In the above equations *C*<sup>1</sup> and *C*<sup>2</sup> represent dry above ground biomass (D-AGB) in ton/ha. Tables 1 and 2 show the cross-correlation/covariance of the factors entering in the first and second stage of regression, respectively.

The D-AGB values are now regressed against the results obtained from the previous stages shown as *C*<sup>1</sup> and *C*<sup>2</sup> bilinear expressions (stage 3). This last regression ends up with the estimation of four parameters *mi*(*i* = 1, ... 4) and the final working formula for D-AGB on a daily basis is given by the following Equation (3) in an implicit form, since *C*<sup>1</sup> and *C*<sup>2</sup> are functions of the attributes LAI, *h*c, and cumETc:

$$\text{D-AGB(ton/ha)} = m\_1 + m\_2 \cdot \text{C}\_2 + m\_3 \cdot \text{C}\_1 + m\_4 \cdot \text{C}\_1 \cdot \text{C}\_2 \tag{3}$$

where *m*<sup>1</sup> = 0.0082, *m*<sup>2</sup> = 1.11, *m*<sup>3</sup> = −0.12, *m*<sup>4</sup> = 0.0032.

The above algorithm is exemplified by a flowchart diagram shown in Figure 1.

**Figure 1.** Flow chart showing the procedure of establishing the new empirical equation for daily estimation of dry above ground biomass (D-AGB).

In Figure 2a, the iso-lines of (D-AGB) derived from Equation (1) as a function of LAI, plant height (*h*c), and D-AGB measurements, through curve interpolation lines, respectively, on a daily basis, are presented. Similarly, Figure 2b shows the results of the second stage D-AGB = *f*(LAI,cumETc). Higher sensitivity showed the LAI-hc correlation rather than the one of LAI-cumETc for the D-AGB factor.

**Figure 2.** The iso-lines of (D-AGB) derived from Equations (1) and (2). (**a**) first stage (**b**) second stage.

Tables 1 and 2 show the cross-correlation/covariance of the factors entering in the first and second stages of regression, respectively.

**Table 1.** Cross-correlation/covariance between LAI, *h*c, and D-AGB from the first stage of regression.


**Table 2.** Cross-correlation/covariance between LAI, cumETc, and D-AGB from the second stage of regression.


It is obvious from the Tables 1 and 2 that the strongest correlation exists between LAI and *h*<sup>c</sup> (0.996) and that all three attributes; LAI, *h*c, and cumETc, are also strongly correlated to D-AGB (all correlation coefficients are above 0.92, see Tables 1 and 2).

The regression equations between daily simulated D-AGB values against the experimental and the cross-correlation coefficient (R2) are shown in Table 3 for the 2015 and 2014 cultivation periods, respectively.

A statistical analysis was further performed in order to provide quantitative indices to our estimates. For this purpose, the following statistical indices were evaluated [19,20]: (i) Mean bias error (MBE), (ii) Variance of the distribution of differences *s*<sup>2</sup> *<sup>d</sup>* which expresses the variability of (*<sup>P</sup>* <sup>−</sup> *<sup>O</sup>*) distribution about MBE, (iii) Root mean square error (RMSE), (iv) Mean absolute error (MAE), (v) Index

of agreement, *d*, [20], where *n* is the number of cases. *O* denotes the experimental values of D-AGB measured during the 2014–2015 cultivation periods for all irrigation treatments (*I*100, *I*75, *I*50, and *I*25). *P* denotes the simulated values as these are estimated by the proposed methodology. All the above mentioned relevant statistical indices are provided in Table 3.


**Table 3.** Summary statistics of daily dry above ground biomass (D-AGB) tested against the reference method.

#### **3. Results and Discussion**

In Figure 3, the development of D-AGB, both measured and simulated during cultivation period 2015 expressed in days after planting (DAP), is presented. As it is depicted in Figure 3a, the simulated and experimental curve interpolated lines almost coincided for the 100% treatment because the model has been calibrated for this treatment and cultivation period. Figure 3b–d depicts the 75%, 50%, and 25% treatments for the 2015 cultivation period, respectively. It is obvious that the predictions by the model for the 75%, 50%, and 25% treatments give results very close to the measurements for the 2015 cultivation period.

**Figure 3.** Relationship between dry above ground biomass (D-AGB), (ton/ha), and days after planting (D-AGB ) for 2015 cultivation period and PR91M10 hybrid. The (**a**), (**b**), (**c**), and (**d**) parts depict the 100%, 75%, 50%, and 25% of ETc water treatments, respectively.

In Figure 4, the development of D-AGB both measured and simulated curve interpolated lines during cultivation period 2014 expressed in days after planting (DAP), are presented. For the 2014 cultivation period, all four figures (Figure 4a–d)were used for verification purposes. Figure 4a–d depicts the 100%, 75%, 50%, and 25% treatments, respectively. From Figure 4a at 100% treatment, it can be assumed that measured and simulated curve interpolated lines were very close, and at DAP 75 the model predicted D-AGB 4.951 ton/ha, while experimental D-AGB for the DAP 75 for the non-water stressed soybean was 4.385 ton/ha. Similarly, the 75%, 50%, and 25% water treatments were perfect fitted until 55 DAP, approximately, for the 2014 cultivation period. However, the response of the plant to the water stress mechanism is a fairly complex process involving both biophysical and biochemical functions that could differentiate predictions of the experimental observations. This induces the differences after DAP 55 for the water stressed treatments in 2014 cultivation period (Figure 4b–d) and for the 2015 water stressed treatments (Figure 3b–d).

**Figure 4.** The relationship between dry above ground biomass (D-AGB), (ton/ha), and days after planting (DAP) for 2014 cultivation period and PR91M10 hybrid. The (**a**), (**b**), (**c**), and (**d**) parts depict the 100%, 75%, 50%, and 25% of ETc water treatments, respectively.

#### **4. Conclusions**

For the first time, an already existing empirical methodology for the prediction of reference evapotranspiration (ET0) used with crop geometrical characteristics (LAI, *h*c) and cumETc as inputs was used in order to predict daily D-AGB. The statistical analysis showed very satisfactory adjustment of the experimental and simulated values, especially for the non-water stressed treatments of the 2014 cultivation period.

Further experimentation for different regions and a wider range of D-AGB values is needed in order to verify the goodness of fit, for the parameters used in the methodology, in different climate regimes, and for more cultivation species. An important advantage of the methodology that has been followed, in addition to the use of three readily measured fundamental parameters (cumETc, LAI, *h*c), is that the model can easily be calibrated (different coefficients) for any crop and in any climatic environment.

However, a more complex algorithm could be set using more environmental attributes of the soil-plant-atmosphere system, which might be adjusted better to the simulated values than the experimental ones, especially for the plants under non-water stress or under irrigation deficit.

**Author Contributions:** Conceptualization (S.A. and I.A.); Investigation (C.V.); Methodology (S.A., I.A., C.V.); Supervision (I.A.); Writing—original draft (C.V.); Writing—review & editing (I.A., S.A., C.V.). All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was supported by the Alexander S. Onassis Foundation.

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

#### **References**


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

## *Article* **Decision-Making Process for Photovoltaic Solar Energy Sector Development using Fuzzy Cognitive Map Technique**

**Konstantinos Papageorgiou 1,\*, Gustavo Carvalho 2, Elpiniki I. Papageorgiou 3,4, Dionysis Bochtis <sup>4</sup> and George Stamoulis <sup>1</sup>**


Received: 25 February 2020; Accepted: 16 March 2020; Published: 19 March 2020

**Abstract:** Photovoltaic Solar Energy (PSE) sector has sparked great interest from governments over the last decade towards diminution of world's dependency to fossil fuels, greenhouse gas emissions reduction and global warming mitigation. Photovoltaic solar energy also holds a significant role in the transition to sustainable energy systems. These systems and their optimal exploitation require an effective supply chain management system, such as design of the network, collection, storage, or transportation of this energy resource, without disregarding a country's certain socio-economic and political conditions. In Brazil, the adoption of photovoltaic solar energy has been motivated not only by the energy matrix diversification but also from the shortages, problems, and barriers that the Brazilian energy sector has faced, lately. However, PSE development is affected by various factors with high uncertainty, such as political, social, economic, and environmental, that include critical operational sustainability issues. Thus, an elaborate modelling of energy management and a well-structured decision support process are needed to enhance the performance efficiency of Brazilian PSE supply chain management. This study focuses on the investigation of certain factors and their influence on the development of the Brazilian PSE with the help of Fuzzy Cognitive Maps. Fuzzy Cognitive Map is an established methodology for scenario analysis and management in diverse domains, inheriting the advancements of fuzzy logic and neural networks. In this context, a semi-quantitative model was designed with the help of various stakeholders from the specific energy domain and three plausible scenarios were conducted in order to support a decision-making process on PSE sector development and the country's economic potential. The outcome of this analysis reveals that the development of the PSE sector in Brazil is mainly affected by economic and political factors.

**Keywords:** Fuzzy Cognitive Maps; photovoltaic solar energy; scenario analysis; decision-support; energy management

#### **1. Introduction**

Fuzzy cognitive maps (FCMs) were introduced by Kosko (1986) as a methodology for modelling complex, dynamic systems [1]. In terms of structure, FCM is a directed weighted graph with feedback loops, which consists of nodes (concepts) and strength connections (relationships) among nodes [2]. An FCM models the system in the form of a graph, introducing causal relations between concepts and

represents the system's behavior using the accumulated knowledge, perceptions, experiences or beliefs from experts, stakeholders or both [3,4], usually expressed in natural language.

Over the last years, FCMs have attained great attention and popularity since they are able to implement modelling, analyzing and simulating tasks, as well as test the influence of parameters and predict the behavior of the examined system [2]. Furthermore, they include learning capabilities and characteristics as they can learn from historical data and previous knowledge to overcome the subjectivity of experts' opinions [5–7]. In other words, it is a semi-quantitative method that uses qualitative knowledge from experts and stakeholders to form FCM-based models and further increase system understanding and reduce uncertainties.

The task of simulation is incorporated by FCMs and can help managers to increase system understanding and reduce uncertainties about how the studied phenomenon may respond to changes in the main system drivers [8].

Among other fields, FCMs are also used for scenario planning even if they were not initially used for this purpose [9]. Kok [10], Van Vliet et al. [11] and Jetter and Schweinfort [12] have introduced FCM-based scenarios to improve the quality of scenarios and increase the robustness of the scenario development step.

Over the last two decades, renewable energy has been considered a vital candidate in global power demand [13], and thus, its exploitation needed to be well explored and studied. However, due to several limitations, such as environmental problems and because of the high uncertainty that these kinds of systems present, the renewable energy sector is considered highly non-linear and complex. In this context, Fuzzy Cognitive Map is a method that can illustrate the uncertainty causality [14] and becomes a powerful tool to model such complex systems that involve many factors [2]. Due to this fact, FCM has been recently deployed in a renewable energy domain with great perspectives.

Taking one step further, FCMs have recently implemented management and decision-making tasks in the Energy sector [9,15]. In particular, in [15], an FCM approach was applied to predict and analyze the consequences of Brexit for the Energy-Water-Food (EWF) demand, examining different Brexit scenarios. FCMs can deal with complex problems, vast domains of inherent uncertainties [16] and processes like decision making that are based on human reasoning process [2,17]. In [18] and [19], special attention has been given to FCMs for management, decision-making and prediction purposes in the energy context. Specifically, the authors in [18] have proposed an algorithm based on Artificial Neural Networks (ANN) to create a model for monthly gasoline consumption estimation. The produced results could help officials and producers for better planning and marketing as well as optimizing gasoline production. An FCM has been proposed in [20] as a powerful inference tool to model the complex system of the electricity market. In this context, the authors describe the structure of FCMs and examine the way they are used to model and simulate competitive trading markets. Thus, certain approaches will be selected to further increase the understanding of how such markets behave in different circumstances. A Fuzzy logic-based framework was investigated in [21] for anomaly detection in energy consumption data where uncertainties inherent to forecasting, have been tackled. Regarding the renewable energy domain, a socio-economic analysis was conducted in [22], in which FCMs were employed to identify potential consequences over the mix of fossil fuel and renewable energy power sectors, as well as how certain policies that regard the energy domain can affect a country's economic growth. We also need to report on [23], in which an FCM methodology was applied for assessing energy efficiency strategies and electricity demand prediction, whereas Nikas et al., dealt with policy-making in the Renewable Energy sector, employing a stakeholder-driven approach based on FCMs [24]. Alipour et al., in a more recent study, explored the FCM capabilities on identifying and analyzing unpredictable factors that have a direct impact on solar photovoltaics in an uncertain socio-economic, political and technological environment [25]. It is worth mentioning that this is the first research study that implements FCM-based scenario analysis in the renewable energy field.

In order to tackle the vast amount of data, various parameters, modelling, learning, and simulation tasks of complex systems, researchers need a simple but powerful software tool that could help them in this direction. As found on the related literature, several tries occurred by researchers for developing a software product with simple or advanced functions that could perform either designing and analyzing tasks or even more complex functions like learning and simulation on FCM-based systems. Reviewing the existing literature, tools like Mental Modeler [26], FCMapper [27], ESQAPE [23], FCM Expert [28], FCM-Analyst [29], JFCM [30], and FCM Designer [31] can offer various capabilities and have their own unique specifications. However, MentalModeler [26] is the only web-based tool for scenario analysis and environmental modeling. This tool provides an opportunity to incorporate different types of knowledge into environmental decision-making, define hypotheses to be tested, and run scenarios to determine perceived outcomes of proposed policies.

However, a gap is noticed in the relevant literature, regarding the existence of a more user-friendly tool that could provide extra capabilities for aggregating FCM models and learning them using supervised and unsupervised algorithms. Considering this gap, the team of Papageorgiou E. has developed a new web-based tool, called FCMwizard [25], that incorporates these functionalities for helping research community to freely experiment with them and exploit all the capabilities offered. Some preliminary work with the use of the FCMWizard tool has been recently made in the environmental domain, dealing with decision making problems such as (i) the economic sustainability of poor women in rural areas [32] and (ii) climate change adaptations [33].

In this work, the FCM technique is applied to perform a decision-making process for PSE sector development, providing at the same time the researchers and the FCM community, wth an efficient simulation tool called FCMWizard, to conduct policy making simulations in their own case studies. In particular, the authors explore relevant FCM-based scenarios for the Brazilian Photovoltaic Solar Energy Sector (PSE) employing the FCMWizard tool. The PSE arises as one of the most promising renewable energy sources to replace fossil fuels [34], since climate and environmental changes, described by the scientific community [35], are increasing continuously. This happens as a result of the increasing energy demand that stems from technological progress and other advancements in human development. Scenario planning has proven to be the appropriate technique for institutions that are immersed in this PSE sector and need to make complex decisions, taking into account various uncertainties and risks, such as socio-political, environmental and technological.

The following tasks emphasize the novelty of this research work:


The structure of the present paper is as follows: Section 2 briefly describes FCMs, while Section 3 presents the problem of PSE development in Brazil and the building process of the model. In Section 4, FCMWizard tool is adequately presented and its decision-making functionalities are reported accordingly. The scenarios development and analysis regarding the Brazilian PSE with the use of the FCMs are provided in Section 5. The discussion of the results is laid in the same section. The last section is devoted to the conclusions of this study.

#### **2. Fuzzy Cognitive Maps**

FCMs are fuzzy structures that combine fuzzy logic and Neural Networks, which incorporate knowledge acquisition and indicate the causality among concepts of complex systems that are characterized by uncertainty. FCMs were first introduced by Kosko [1] and since then they have been used to represent human knowledge and reasoning, while they serve as a method for modelling dynamic decision systems. Furthermore, learning techniques were later included to train FCMs and identify appropriate weights for its interconnections.

FCMs as we know them, are fuzzy signed directed graphs with feedback. They consist of concepts that are represented in the graph by nodes. Concepts are interconnected with directed arcs that are labelled with fuzzy values in the interval [0, 1] or [−1*,* +1] and show the strength of impact between the concepts [36]. Typically, an FCM with *n* nodes can be represented in the form of an adjacency matrix W*nxn*, where each element *wij* of the matrix takes values within [−1, +1].

The causal relationships (weights) between concepts have three possible types.


The FCM converges to an equilibrium point using one of the following calculation rules, (1)–(3), when values X*<sup>i</sup>* have been assigned to the concepts *Ci* [37]:

$$\text{Kosko}: X\_i^{(\kappa+1)} = f \left( \sum\_{\substack{j=1 \\ j \neq i}}^n w\_{ji} \times X\_j^{\kappa} \right) \tag{1}$$

$$\text{Modified Kosko}: X\_i^{(\mathbf{x}+1)} = f \left| X\_i(t) + \sum\_{\substack{j=1 \\ j \neq i}}^n X\_j(t) \cdot w\_{ji} \right| \tag{2}$$

$$\text{Rescale}: \ X\_i^{(\mathbf{x}+1)} = f\left(\mathbf{2} \times X\_i^{\mathbf{x}} - 1\right) + \sum\_{\substack{j=1,\ j \neq i}}^n w\_{ji} \times (2 \times X\_j^{\mathbf{x}} - 1)\tag{3}$$

where *<sup>X</sup>*(κ+1) *<sup>i</sup>* is the value of concept *Ci* at simulation step <sup>κ</sup> <sup>+</sup> 1, *<sup>X</sup>*<sup>κ</sup> *<sup>j</sup>* is the value of concept *Cj* at the simulation step κ, *wij* is the weight of the interconnection between concept *Ci* and concept *Cj*, κ (Greek letter Kappa) is the interaction index and every simulation step, and f(·) is the threshold (activation) function. The function is selected to retain the values within the range [0, 1] or [−1, 1]. Equation (2) is the original Kosko activation function (1) that takes past values (or history) of a concept into account [38]. Generally, there are four most used transformation functions (4)–(7): (a) bivalent, (b) trivalent, (c) sigmoid, and (d) hyperbolic tangent.

$$\text{Bivalent}: f(x) = \begin{cases} 1, & x < 0 \\ 0, & x \ge 0 \end{cases} \tag{4}$$

$$\begin{array}{c} \text{Trivial:} \quad f(\mathbf{x}) = \begin{cases} \ 1, \mathbf{x} > 0 \\ \ 0, \mathbf{x} = 0 \\ -1, \mathbf{x} < 0 \end{cases} \end{array} \tag{5}$$

$$\text{Signpoint (linear)}: f(\mathbf{x}) = \frac{1}{1 + e^{-\lambda \mathbf{x}}} \tag{6}$$

$$\text{Hyperbolic tan get}: f(x) = \tanh(\lambda x) \tag{7}$$

where λ is a real positive number (λ > 0), which determines the steepness of the continuous function f and x is the value *X*<sup>κ</sup> *<sup>i</sup>* on the equilibrium point. At each step, the value *Xi* of a concept is influenced by the values of concepts connected to it and it is updated according to the inference rule. This process continues until the system converges, which indicates that the difference between two subsequent value of the outputs values must be equal or lower to ε (epsilon, ε = 0.001).

#### **3. PSE in Brazil**

#### *3.1. Problem Statement*

Environmental degradation caused by deforestation, gas emission, and environmental pollution has caused countries to rethink their electric systems. As a solution to the problems caused by the unsustainable exploration of fossil fuels, renewable energies have become the focus of attention of a broad range of agents. However, in Brazil, unlike other countries, the investment in photovoltaic solar energy was driven by different reasons, such the increase in the cost of energy produced in thermoelectric plants and the emission of greenhouse gases by burning of fossil fuels [39,40].

Among other solutions, PSE shows up as a viable and necessary source of electrical energy. Given that it is located mostly in the intertropical region and has excellent potential for solar energy utilization throughout the year, the Brazilian solar radiation is higher than the European for almost its entire territory [41]. Based on Renewable Energy Country Attractiveness Index [42], Brazil is eighth in the raking of attractiveness when speaking of photovoltaic solar energy. In 2014, Brazil had its total installed to power up to 15 MW, and in 2015, the same number surpassed 32 MW. Statistics of the Mines and Energy Ministry [43] indicated that by 2018 Brazil should be between the 20 countries with the most prominent generation of solar energy.

On the other hand, the high cost of acquisition of photovoltaic and their low conversion efficiency from solar irradiation to electrical energy are the main impediments to the large-scale diffusion of these systems [44,45]. As a solution, the government can understand the dynamic of PSE development from other countries and propose policies that improve the use of solar energy in an urban environment [46]. Given that the future of the PES depends on several issues, it needs to be planned and controlled. For this purpose, scenario planning is a fundamental tool [47].

#### *3.2. Model Development*

In order to provide a quasi-quantitative model for scenario planning, regarding the studied problem of PSE development in Brazil, the FCM development process was conducted in six steps. Figure 1 offers a visual representation of these steps for better understanding of the procedure and for the convenience of the readers.

**Figure 1.** The steps of FCM development process.

#### Step 1: Obtaining concepts from participants.

The first step deals with the determination of the most essential independent variables that define the examined problem and affect the dependent variable. Participants are asked to determine the concepts of the system but only the most important in order to avoid designing a system with a vast number of variables that would be difficult to study.

Step 2: Qualitative assignment of cause-effect relationships between concepts.

In the second step, participants are asked to assign values on the scale of 1–10 and to determine whether there is a positive or negative cause-effect relationship for the weighted interconnections among the depended and most significant variables that constitute the FCM model. Ten (10) denotes the highest strength and one (1) the lowest. The sign + (plus) denotes a positive influence, whereas the sign – (minus) denotes a negative influence that one node has to another.

Step 3: Weights Normalization (coding into an adjacency matrix).

Weights given to each link were then normalized between 0 and 1 (as described in Table 1), considering positive and negative values for coding into the adjacency matrix [48]. So, we substitute the qualitative values assigned by the stakeholders for expressing the degree of influence with the corresponding quantitative values, in order to define the weight matrix of each stakeholder, as presented in Table 1. For example, the qualitative value 10 (that denotes the highest strength) is substituted by the quantitative weight value 1.

**Table 1.** Normalizing Qualitative (linguistic) values to Quantitative values for the weighted interconnections.


Step 4: Producing individual FCMs from each group.

In this step, every group of the participants was asked to construct an individual FCM by defining the primary variables and determining the weights values of all causal relationships. This process also includes the identification of the decision concept, which is a dependent variable for the problem under investigation. As we aim to analyze the Brazilian Solar Photovoltaic Energy Sector, the dependent variable was set to be the "Development of PSE in Brazil". The procedure that was followed, is described below in detail, as the authors want to highlight the importance of all the steps taken for the success of this study.

Step 5: Groups Aggregation producing an Overall FCM.

In this step, all individually coded cognitive maps from the three groups were aggregated and an overall group FCM (Collective-FCM) was produced that includes all the concepts from all individual cognitive maps. This collective FCM represents the perception of all the stakeholders and is enriched with the knowledge of all stakeholders involved.

Step 6: Visualization of collective FCM.

After the aggregation process, that was based on the weighted average method, the Collective-FCM was analyzed using the FCMWizard software tool (version 1.0, E.I. Papageorgiou, Larissa, Greece). Since the tool includes modelling and visualization capabilities, a visual representation of the condensed FCM model was created by FCMWizard, which specifically illustrates the concepts and all the connections between them. The graphical representation of the collective FCM is presented in Section 4.1, where an overview of the FCMWizard tool is presented.

Following the steps above, the authors carried out the project with the help of researchers of the University of São Paulo (USP) and University of Brasília (UNB) specialized in the solar photovoltaic energy. They along with other specialized stakeholders (such as specialists from the National Electric Energy Agency (government), the Brazilian Solar Energy Association (professionals)) were the primary sources of data and information for the development of this research.

Interviews were conducted individually or with pairs of specialists from the National Electric Energy Agency (government), the World Wildlife Fund (NGO), Brazilian Solar Energy Association (professionals), and researchers from the University of São Paulo and Brasília (specialists). A workshop was also carried out with a group of eight graduate students of the Institute of Energy and Environment of the University of São Paulo, IEE/USP, to consolidate the FCMs.

Specifically, a pretest was first carried out with a potential consumer with a business background, through an individual interview, during which the dynamics of FCMs and the proposed method were explained. The respondent had trouble defining the primary variables (concepts) and establishing the exact weights (from −1 to 1) of the causal relations. After proper guidance and following the contents of Table 1, in the end, he managed to develop a potential FCM.

The second FCM was constructed after interviewing a specialist in SPE from the University of Brasília (UNB). He responded with great enthusiasm and a consensus was reached for 12 interrelated concepts, establishing at the same time the causal relations and respective weights.

In the next FCM, the proposed method was investigated with the help of a group of SPE specialists from IEE/USP. In a workshop with eight participants, the dynamics of FCMs and the proposed method were explained. These specialists identified 13 concepts and established their connections and weights. In this type of approach (workshop), the most prominent advantage came from the debate among the participants. Instead of creating individual FCMs, the authors suggested a collaborative process for scenario planning where the participants cooperated to produce a more productive and accurate knowledge in the form of a consolidated FCM. The weight value for each interconnection is calculated as the average of all values that each participant gave for the corresponding interconnection, taking into consideration Table 1. However, the FCM achieved by the group of specialists from USP, presented in Figure 2, was not a significant improvement regarding complexity and robustness in comparison with the others.

**Figure 2.** FCM from a group of the specialist of IEE/USP (average weight values).

The FCMs elaborated from the information of the specialists, and stakeholders were integrated into a single FCM. Table 2 presents the list of the concepts contained in this collective FCM, with a short description of them. They identified 29 concepts in total and established the connections and weights among these concepts. Also, they determined the 10 most central concepts (concepts in bold in Table 2) in terms of significance in the decision making process through scenario analysis. Moreover, the weights of the connections were equal to the average of the weights established by the specialists and stakeholders.


#### **Table 2.** List of concepts integrating the collective FCM.

In this context, a 29-nodes FCM model was produced from the experts/stakeholders and its adjacency matrix (see Figure 3) was created accordingly, including the causal relationships between the concepts and their weight.


**Figure 3.** Adjacency matrix.

#### **4. Overview of the FCMwizard Tool**

This section provides an overview of the FCMWizard tool in which the functionalities of this software are presented. In order to illustrate the capabilities of this tool on modelling and simulating a dynamic and complex system, the authors use the case study of the renewable energy sector of Brazil. Among the functionalities that site on the main menu (see Figure 4), there are three main modes for constructing an FCM, namely: (i) Expert mode wherein experts and/or stakeholders can construct manually a FCM for a real problem; (ii) Data-based Learning mode wherein the FCM model is constructed automatically using given data; and (iii) Merge mode, where different individual FCMs can be combined to provide an augmented and collective FCM model in order to generate aggregated system complexities.

**Figure 4.** FCMWizard Main menu, with FCM construction and FCM Learning modes.

FCMWizard offers a Main menu with 2 main functionalities: FCM construction and FCM Learning. FCM Learning includes three learning modes (see Figure 4): (i) Hebbian Learning (NHL): weights are

fine-tuned by data assuming that all the concepts are synchronously triggered at each iteration and they synchronously change their values. Constraints on concepts and weights can be considered for updating weighted connections. (ii) Data-driven Hebbian Learning (DD-NHL): an extended NHL method through historical data for the input concepts to fine-tune weights, and (iii) Differential Hebbian Learning (D-NHL): This algorithm correlates the changes in the concepts to modify their weights [2].

Next, the authors will illustrate how the FCM Wizard tool can be used to construct an FCM model and visualize its results for the given problem of PSE in Brazil. Moreover, the structural analysis, the Hebbian-based Learning and the scenario analysis will be further described, along with relevant figures and tables, produced by the tool.

#### *4.1. Expert-Based Construction*/*Modelling*

Without a deeper knowledge of the mathematical foundations of the methodology and regardless of the type of application domain, users can create new concepts and causal relationships between them, using influence weights, to design a new model. The "Expert Mode" under the main menu "Construct FCM", provides users with the ability to easily design an FCM model by adding new nodes when the "+" button in the "Model Design" tab is selected. Then, new directed lines can be drawn representing the causal relationships between nodes, in which a green direct relationship denotes a positive weight and a red inverse relationship denotes a negative weight.

So, the FCM model that was designed by the stakeholders can be easily constructed, as illustrated in Figure 5, using the FCMWizard tool, following the process described above. This consensus FCM model comprises the concepts in the form of nodes (C1, C2, ... , C29) and directed lines (green or red) with the corresponding weight (positive or negative). Moreover, this FCM model can be directly designed by the FCMWizard, after properly importing the adjacency matrix (see Figure 3) in the tool, as provided by the Stakeholders. The weight matrix, as illustrated in the tab "weight matrix" of the tool, is presented in Figure 6. Each number of this weight matrix represents the causal relationship between the concepts that are sited in the corresponding row and column of this matrix.

**Figure 5.** Consensus FCM model for PSE in Brazil designed in FCMwizard tool.


**Figure 6.** Screenshot of weight matrix of the final consensus FCM model produced in FCMWizard.

#### *4.2. Structural Analysis*

This software module encapsulates the ability to calculate various graph theory indices like the total number of concepts, number of connections, connection-to-concept ratio, outdegree, indegree, type of variables (receiver, transmitter or ordinary), degree centrality, betweenness centrality, closeness centrality, complexity ratio, density, and hierarchy index [33,49]. In Figure 7, a snapshot of the first eight concepts' graph indices is illustrated. Centrality is considered a significant index, so that a concept with a high degree of centrality has an important role in the cognitive map [8,49]. Centrality is calculated by the sum of the corresponding absolute indegree and outdegree causal weights [50]. Table 3 gathers the two most important indices of centrality (Degree and Betweenness) for all 29 concepts, as calculated by the tool. As observed from Figure 7 and Table 3, the concepts C1, C4, C6, C8, C11, C14, C16, C18, C24, and C25 present a higher degree of centrality than the rest of the concepts, so they can have a significant role on the examined energy system.

Complexity, density and hierarchy index too, are among significant structural indices for analyzing the graphical structure of FCMs.

From the constructed FCM model (Figure 5), it emerges that concepts with the most substantial number of connections are: Public knowledge on SPE (C4), Public opinion (C5), Energy demand (C8), Government incentives (C11), Private sector involvement (C12), Environmental pollution (C17), Energy price (C18), and Economic situation (C25). These key concepts are also directly and strongly connected to the objective concept C1 (development of solar photovoltaic energy). Hence, in terms of their degree of centrality and connectedness (see Figure 7 and Table 3), concepts C4, C5, C11, C17, and C18 were chosen by the researchers to become the base of the scenarios since they can strongly affect the behavior of the system.


**Figure 7.** Screenshot of the "Structural Analysis" tab.



#### *4.3. Scenario Analysis*

FCMWizard also provides users with the ability to perform "what-if" scenario analysis after they specified the desired parameters of the FCM model to simulate [51,52]. The initial stimulus state vector, the inference rule's type, the transfer function, its learning parameter, and the number of iterations or the convergence step need to be properly defined. Users are also offered the open/closed lock option to set and keep the value of a concept unchained throughout iterations, as being "clamped" [25].

An example of a simulation model run by the tool for the examined case is presented below and regards the investigation of the impact that concept C11-"Government incentives" has on other concepts as well as on the examined system, overall.

Before the simulation process takes place, a baseline scenario needs to be conducted, so in the next step there will be a comparison between the results of the three scenarios and those resulting from the state of the model where all the initial values of concepts are zero (baseline scenario). In Figure 8, an exemplary table with the first 20 concepts values in all iterations steps, regarding the baseline scenario, is presented, while Figure 9 depicts the corresponding graph showing the steady state of every concept comprising this FCM model.


**Figure 8.** Screenshot of the FCM iterations steps in the baseline scenario from FCMWizard.

**Figure 9.** Screenshot of the graph regarding concepts' values in the Baseline scenario.

Let's now consider that concept "Government Incentives" (C11) is activated for the considered PSE problem, and the process of how this concept affects the other concepts needs to be investigated. For performing the scenario analysis, authors need to define a number of parameters like the inference rule, the transformation function, Lambda parameter, and either the number of iterations or the convergence step. All these simulation options and learning parameters are offered by the tool. In particular, it supports three FCM inference rules, four threshold functions, besides the possibility of customizing other parameters too [8].

In Figure 10, a screenshot of the options provided by "FCMWizard" is illustrated, regarding the parameters' configuration for the simulation process. Specifically, the initial value of the concept C11 was clamped to 1, the Modified Kosko's activation rule and the Sigmoid transformation function were selected, as well as the Lambda parameter and the convergence step were set to 1 and 0.0001 respectively.

The produced values for every concept of this scenario analysis example in the form of table and graph are depicted in the following figures (see Table 4 and Figure 11). The rest of the results are presented in Section 5.2, where the variety of functionalities and capabilities of the FCMWizard tool is further revealed.

#### *Energies* **2020**, *13*, 1427


**Figure 10.** Screenshot of the Scenario analysis mode by FCM Wizard.



**Figure 11.** Graphical plot of convergence performed by FCMWizard.

Through this process, FCMWizard is proven to be a useful tool for designing FCM models and building simulations for different scenarios. It provides easy-to-use but powerful capabilities in modelling, analyzing and simulating complex systems with high uncertainty, and thus, it can be applied by policymakers for strategic decisions in various scientific domains. It should be noted that this web-based software is freely available to users upon a password request and is under continuous development.

#### **5. Applying FCMWizard for Scenario Analysis**

For a deeper understanding of the examined problem, such as the concept's behavior, their relations, and the extent to which one concept has an impact on others, an FCM simulation process is needed. This study employs the FCMWizard tool to simulate and analyze the problem [25]. This process entails the multiplication of the input vector with the adjacency matrix, and in the output vector a squashing function (Equation (2)) is applied each time as a threshold function, until the system reaches a stable state.

The process of simulation is performed by "clamping" the initial value of certain key concepts each time. This outcome is compared against a baseline scenario where the system reaches the steady state. After analyzing the deviations of concepts values between the baseline steady state and the outcome of the simulation process, researchers can interpret the impact of the key concepts on the system on a quantitative basis.

#### *5.1. Scenario Development*

The first established approach in scenario planning was the selection of the most important concepts (called decision concepts). Among the concepts selected to construct the studied FCM model (see Table 2), the researchers identified five concepts that could well affect the behavior of the system. These concepts (C4, C5, C11, C17, and C18), given by the programme participants and implementers, were selected as they were among the concepts with the highest centrality, having both in/out-degree values, while they are strongly connected with the objective of this scenario planning, which is the development of PSE (concept C1).

For the purposes of this study, three scenarios were developed through plausible combinations of the aforementioned decision concepts. To check the behavior of the system as a whole, we tested a set of combinations in which the concepts' values were activated and clamped to one. Thus, we create three scenarios: "Instability Scenario", "Disastrous" and "Government Incentives".

The selected scenarios with their concepts are briefly presented in the following table (Table 5).


**Table 5.** The key concepts of each scenario.


The first scenario, namely "Instability Scenario", reflects the recent history of Brazil, in which the country is undergoing economic and political instability. The current economic crisis in Brazil began in mid-2014, and this economic crisis was accompanied and intensified by a political crisis [53]. As a result of these crises, Dilma Rousseff, president of the time, suffered impeachment in 2016 [54]. The crisis also generated unemployment, which reached its peak in 2017 [55]. In response, the price of energy has risen dramatically. Consecutively, the concept "Energy price" was chosen to be the only concept examined in this scenario and so, its value is clamped to 1.

The "Disastrous" scenario reflects the shortage of electricity in Brazil, the blackouts, as well as some other environmental disasters, namely Mariana and Brumadinho [56]. Brazil has had to activate thermoelectric plants regularly, but this backup capacity was insufficient, and rationing had to be imposed. To avoid another "blackout crisis," the government decided to create a set of initiatives to encourage the growth of PSE, through the establishment of subsidized credit lines and tax exemptions. In 2015 and 2019, the mining dams named Brumadinho and Mariana, collapsed. Part of the demand for energy migrated to the PSE and the price of electricity provided by the government decreased. This scenario examines the impact the concepts "Public knowledge on PSE" and "Government incentives" have on the system dynamics and were accordingly activated and clamped to 1.

In the third, "Government Incentives" scenario, Government offers financial assistance to private businesses making investments through the use of economic incentives. Incentives can include tax abatements, tax revenue sharing, grants, infrastructure assistance, no or low-interest financing, free land, tax credits, and other financial resources. Thus, there seems to be more private sector involvement leading to a greater market competition, resulting in an overall energy price decrease. This could make the public feel more comfortable to consume more energy. In this scenario, the concept "Government Incentives" was clamped to 1, so the researchers can estimate through the simulation process the degree this concept affects the examined energy system.

#### *5.2. Scenario Analysis for PV Solar Energy Sector*

As it was reported in the "Scenario Analysis" section, the simulation process started after the baseline scenario was performed. Then, scenario analysis performed simulations for the selected three scenarios (Table 5). In this context, scenario 1 (S1) was devoted to increasing the decision concept C18—"Energy price" by "clamping" it to one. Scenario 2 (S2) studied the effects of the decision concepts C4—"Public Knowledge on PSE", C5—"Public opinion", C11—"Government Incentives" and C17—"Environmental pollution" by clamping the values of these concepts to one. In the end, Scenario 3 (S3) investigated how the increase of decision concept C11—"Government Incentives" affects the examined system, when "clamping" its values to one. The results for scenario analysis are illustrated in the following figures. The results are interpreted by comparing the relative change between this baseline scenario and the new steady state.

Figures 12–14 illustrate separately the outcomes of the three scenarios performed (S1, S2 and S3 respectively), with respect to the percentage of change (deviations from the steady state) for each concept. Thus, the authors are able to analyze and further interpret the results produced for the examined renewable energy system scenario, that will help policy makers select the right strategies in this direction. Negative values mean that the concepts have negative acceleration, whereas positive values imply that the concepts have positive acceleration. Worth mentioning is that negative change in the value of development of solar energy (C1) is regarded as a deceleration in deployment and diffusion of solar energy compared to the current state.

The end goal of the examined case study was set to be the development of the solar photovoltaic energy sector in Brazil, which is represented by the concept C1. Therefore, the scenario analysis mainly focuses on the impact that the decision concepts have on this output concept C1—"Development of solar photovoltaic energy", which is considered by the researchers as the main objective of this project. Furthermore, there is an extended discussion about how the rest of the decision concepts involved in scenario analysis, directly affect in a positive or negative way other significant concepts of the examined PSE sector.

#### *5.3. Discussion of Scenario Analysis Results*

The problem of Brazilian PSE sector development is tackled with the technique of FCMs using FCMWizard as a powerful simulation environment for policy making in the renewable energy domain. The researchers have the ability to form certain scenarios with the help of energy specialists, and stakeholders form the specific domain and conduct simulations in order to observe how the examined system responds to changes of certain concepts' values. Casting a thorough look on the results presented in the figures above (Figures 12–14), certain conclusions can emerge for each scenario conducted.

**Figure 12.** Percentage of change (deviation) for all key concepts when concept C18 is activated (S1).

**Figure 13.** Percentage of change (deviation) for all key concepts when concepts C4, C5, C11, and C17 are activated (S2).

**Figure 14.** Percentage of change (deviation) for all key concepts when concept C11 is activated (S3).

The overall observations that were drawn from the figures above, focus on the following points:


Focusing on the first scenario (S1), outcomes showed that an increase in energy price led people to seek new energy sources alternatives and thus, energy demand was decreased, while public awareness of the SPE was increased, given that the payback period decreased. Furthermore, since the companies might have identified a very turbulent and disadvantageous market to exploit, a reduction in the private sector involvement is observed. Without the participation of the private sector, the purchase cost increased.

The conclusions for the second scenario (S2) are briefly presented as follows. Since environmental pollution has sharpened people's ecological awareness and knowledge on PSE, the Brazilian government creates various initiatives and programs to foster the growth of this sector, reducing the payback time of the investment. Part of the demand for energy migrated to the PSE and the price of electricity provided by the government, decreased. The purchase cost decreased accordingly.

In response to the third scenario (S3), where government incentives are considered, economic incentives and other financial resources, given by the government, can possibly cause a significant increase in private sector involvement and greater competition in the PSE market, resulting in a bigger employment capacity. A decrease in energy price is possible and consequently, the country's own inhabitants will feel more confident to consume more energy in their daily life.

Taking into consideration the key concepts of the three conducted scenarios, it emerges that government incentives are mostly able to change the stability of Brazil's Renewable Energy system. After a critical overview of the overall results, it is unlikely the PSE sector develops much, unless there are noteworthy political and economic changes in the status of the country. Moreover, the FCM-based simulation process has shown that Brazil's PSE development is dependent on economic parameters such as the incentives given by the government and is influenced by uncertain factors like public ecological awareness. In addition, the penetration of the private sector in the PSE can affect the dynamics of the system as well as have an impact on the everyday life of Brazilian citizens.

The conducted scenarios highlighted the system's uncertainty and showed to the policy makers how they can develop effective strategies for creating a more robust system that is free from possible future instabilities. The tool settings that are offered, can help users to conduct further scenarios and try different configurations by activating different variables, individually or in a set, as well as by changing the weighted interconnections among them, in order to see how the model reacts in different circumstances.

#### **6. Conclusions**

This study deals with the contribution of the FCM method with the help of an efficient simulation tool called FCMWizard, which has been developed to implement/deploy certain tasks such as modeling experts' knowledge, learning and simulating FCM-based complex systems, in the participatory domain. The main contribution of this work lies in the fact that the authors provide experts and stakeholders that belong to various scientific areas, a free under request, user-friendly and flexible, as regards the number of variables and weights, tool, which can be used anytime and in any discipline or domain, without time or memory restrictions. This software that uses certain learning algorithms, offers to users the ability to easily and intuitively design a model for a particular problem, as well as perform scenario planning for making strategic decisions and policy making in a simple, user-friendly environment, without any computational cost.

The innovation in this work is based on the fact that, decision-makers can apply the Fuzzy Cognitive Map method along with the proposed software tool to model experts' or stakeholders' perceptions and further conduct different scenarios implemented by the FCMWizard tool. In this context, policy-making can be performed in various domains and with different configurations of the tool, regarding the type of inference rules, the number of variables, and the weights values or concepts' selection in the scenario planning, thus, highlighting the generic applicability and usefulness of the new tool.

To show the functionalities and examine the performance of this flexible and interactive web-based software, the authors explored a case study concerning the development of Brazil's Photovoltaic Solar energy sector and a relevant scenario analysis was performed. The main outcomes of this research will help policymakers that belong to the Renewable Energy domain, determine certain strategies for better Supply Chain Management, contributing to optimal exploitation of the PSE sector and Brazil's overall socio-economic growth. In this context, a variety of simulation options and different learning parameters were used for proper configuration of the examined system, in terms of modeling, analyzing and simulating the FCM model. In particular, three FCM inference rules and four threshold functions can be selected by any type of user (experts, stakeholders, policy makers or others), while a number of other useful parameters are offered for customizing simulations.

Certain limitations of this work need to be considered in this section. For example, the provided tool is under continuous development to further include more and advanced functionalities, such as the straightforward use of fuzzy sets for FCM design, interval fuzzy values for weights, other efficient and state-of-the-art fuzzy methods (intuitionistic, type 2-fuzzy sets) for FCM model design, as well as other constraints on the concepts that constitute the aggregated FCM model.

Future work is directed towards further investigation of more scenarios and for different scientific domains. Moreover, for validation purposes, this method could be applied in a respective Renewable Energy domain, where more variables will be considered, while the producing results would be further analyzed and compared to those of the current study. Also, future work is devoted to the development of the data-based construction of FCMs where no qualitative or expert/stakeholder knowledge is available to design an FCM model for scenario analysis.

**Author Contributions:** Conceptualization, K.P.; methodology, K.P. and G.C.; software, K.P.; validation, K.P., G.C., E.I.P. and D.B.; formal analysis, K.P. and G.C.; investigation, K.P. and E.I.P.; resources, G.C.; data curation, G.C. and K.P.; writing—original draft preparation, K.P.; writing—review and editing, E.I.P. and D.B.; visualization, K.P.; supervision, E.I.P., and G.S.; project administration, E.I.P. and D.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

#### **References**


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

## *Article* **A Cloud-Based In-Field Fleet Coordination System for Multiple Operations**

**Caicong Wu 1,2, Zhibo Chen 1, Dongxu Wang 1, Bingbing Song 1, Yajie Liang 1, Lili Yang 1,\* and Dionysis D. Bochtis 3,\***


Received: 16 November 2019; Accepted: 5 February 2020; Published: 11 February 2020

**Abstract:** In large-scale arable farming, multiple sequential operations involving multiple machines must be carried out simultaneously due to restrictions of short time windows. However, the coordination and planning of multiple sequential operations is a nontrivial task for farmers, since each operation may have its own set of operational features, e.g., operating width and turning radius. Taking the two sequential operations—hoeing cultivation and seeding—as an example, the seeder has double the width of the hoeing cultivator, and the seeder must remain idle while waiting for the hoeing cultivator to finish two rows before it can commence its seeding operation. A flow-shop working mode can coordinate multiple machines in multiple operations within a field when different operations have different implement widths. To this end, an auto-steering-based collaborative operating system for fleet management (FMCOS) was developed to realize an in-field flow-shop working mode, which is often adopted by the scaled agricultural machinery cooperatives. This paper proposes the structure and composition of the FMCOS, the method of operating strip segmenting, and a new algorithm for strip state updating between successive field operations under an optimal strategy for waiting time conditioning between sequential operations. A simulation model was developed to verify the state-updating algorithm. Then, the prototype system of FMCOS was combined with auto-steering systems on tractors, and the collaborative operating system for the server was integrated. Three field experiments of one operation, two operations, and three operations were carried out to verify the functionality and performance of FMCOS. The results of the experiment showed that the FMCOS could coordinate in-field fleet operations while improving both the job quality and the efficiency of fleet management by adopting the flow-shop working mode.

**Keywords:** agricultural machinery; fleet management; auto-steering system; collaborative operating system; flow-shop; simulation; field experiment

#### **1. Introduction**

Agricultural fleet management is an important research domain that concerns resource allocation, operation scheduling, routing, communication, and real-time data acquisition [1]. Seamless communications and collaboration across spatial, temporal, operational, and machinery domains are a particular direction toward fleet management in scaled agricultural machinery service organization [2]. Sørensen and Bochtis [3] proposed that the total supply chain of large-scale harvesting has four levels of collaboration, i.e., in-field logistics, inter-field logistics, inter-sector logistics, and inter-regional logistics. In terms of machinery collaboration, Wang et al. [2] proposed that in-field logistics can be subdivided

into three levels, i.e., multi-operation coordination, multi-machinery coordination, and multi-factor coordination [4]. Typical cases of multi-operation collaboration include the master–slave approach [5–7] and leader–follower approach [8–10], which illustrate multiple machines working in parallel for different but related operating tasks. These systems refer to identical operations. Other approaches have been proposed for sequential operations. As an example, in [11], a multi-operation collaboration system was proposed for six sequential operations in a potato production system—namely, bed forming, stone separation, planting, reloading unit, spraying, and harvesting. In parallel, algorithms for optimizing in-field route planning have been developed in various studies [12–15]. Furthermore, the flow-shop working mode of agricultural machinery operation is somewhat more complex [16–19] when multiple units of machinery and multiple operations exist within one field during the same period [20].

For the implementation of all the approaches mentioned above, Real-Time Kinematic–Global Navigation Satellite System (RTK-GNSS)-based automated steering systems can provide automatic guidance for single agricultural machines with one-inch repeatability [21,22]. Auto-steering technology can reduce overlaps and skips and can reduce the dependence upon skilled operators [23]. Therefore, fusion of the efficient flow-shop working mode and advanced auto-steering technologies can improve the field operation efficiency of agricultural machinery. To realize this fusion, the guidance lines followed by the auto-steering systems and the operating progress should be shared among the agricultural machinery units of the fleet in real time, so as to allow the operators to choose the appropriate operating strips for working. It is worth mentioning that different agricultural machinery operations have differing operating strips (in terms of number and size) and, therefore, different guidance lines, since different agricultural machines probably carry different implements with different operating widths. Most auto-steering systems can share guidance lines between different agricultural machinery units for the same operation through U disk copying, radio transfer, and/or server exchange. For instance, Cao et al. [24] built a Web-GIS (Geographic Information System) to realize collaboration among multiple agricultural machinery units, and the John Deere company utilizes an MTG (Modular Telematics Gateway) intelligent gateway and JDLink software to share guidance lines between different agricultural machinery units [25]. To support multi-operation collaboration, Wang et al. [2] designed a cloud-terminal-based coordinative operating system for in-field fleet management. The authors gave the definition of an operating strip to replace the path [26], made use of intervals for strip identification, and proposed an algorithm for strip state updating between adjacent machinery operations. The research demonstrated how to calculate the waiting time or flow step accurately to reduce the delay between adjacent machinery operations [2,19]. However, the validation of the system took place through simulation, while the strip state-updating algorithm provided a feasible—yet not optimized—procedure.

In order to optimize the state-updating algorithm and to validate the functionality of the work presented in Wang et al. [2] under real field experimental conditions, a new collaborative operating system for fleet management (FMCOS) was redesigned and developed, which can support all the operating stages of a field crop and incorporates an upgraded algorithmic process for strip state updating that simplifies the calculation process. The algorithm was tested through simulations using Matlab (MathWorks®) before the development of FMCOS. Finally, three field experiments were carried out to validate FMCOS and to test its functionality and performance.

#### **2. Description of the Fleet Management System**

#### *2.1. System Composition*

Typically, a complete cropping cycle consists of water management, land cultivation, seeding, nutrient management, pest management, and harvesting, as illustrated in Figure 1 [27]. For scaled agricultural machinery cooperatives, adoption of the flow-shop working mode in the phase of land cultivation is needed in order to improve the operating efficiency in a single field. Therefore, FMCOS should not only service the whole production cycle of the field crop, but also support the in-field

flow-shop working mode for some specific operating phases. For this purpose, the auto-steering systems, a fixed GNSS base station, and a database that can manage the parameters of the field area and implements were deployed as a system. The baseline in each field can be used to generate guidance lines for all related machinery based on the width of the implement. FMCOS mainly includes two components: the GNSS-based auto-steering systems installed on the agricultural machinery and the cloud-based operation coordination system.

**Figure 1.** Collaborative application for the production cycle of field crop planting and management.

#### *2.2. Operating Strip Segmentation*

As shown in Figure 2, the field is segmented into strips for each machine. The width of the strip equals the width of the implement. Here, *S<sup>k</sup> <sup>i</sup>* refers to the *<sup>i</sup>*th strip of the *<sup>k</sup>*th operation, *Ok*. Strips are numbered from left to right. *U<sup>k</sup> <sup>j</sup>* are used to describe the machinery units.

A state shows the process progress in a strip by a machine and has four classes: *unready* (0), *ready* (1), *doing* (2), and *done* (3) (Table 1). When the task begins, the state of the strips that belong to the first operation is *ready*, and the machine arranged for the first operation can enter each of the strips for the first operation.

While the fleet is operating in the field, the strip states are updated in real time from the first operation to the last operation in turn. For a single operation, the strip states are updated from the left to the right. When the server gets the updated information of the operating strip from the Android application, the algorithm successively updates the strip state of the current operation according to the prior operating strip state.

**Figure 2.** Expressions of operation tasks, operating strips, and machinery units.


**Table 1.** Definitions of operating strip states.

Figure 3 shows two operations, *Ok*<sup>+</sup><sup>1</sup> and *Ok*. *Ok*<sup>+</sup><sup>1</sup> is the operation after *Ok*, and the state of *Sk*+<sup>1</sup> 3 of *Ok*+<sup>1</sup> is determined by three relevant strips of *Ok*, i.e., *Sk* <sup>4</sup>, *<sup>S</sup><sup>k</sup>* <sup>5</sup>, and *<sup>S</sup><sup>k</sup>* <sup>6</sup>. This is because at the same location in the field, only when the related strips of *Ok* have been finished can the machinery for *Ok*+<sup>1</sup> enter the location. Therefore, the state (*ready* or *unready*) of *Ok*<sup>+</sup><sup>1</sup> is determined by that of *Ok*. In other words, if, and only if, the state of the three strips becomes *done*, the state of *Sk*<sup>+</sup><sup>1</sup> <sup>3</sup> can be updated from *unready* to *ready*.

**Figure 3.** Operating strip state updating.

In general, to update the current strip state, the starting ID (*xk s*) and ending ID (*xk <sup>e</sup>*) of the prior operating strips should be determined. The two IDs and the total number (*n*) of related strips can be calculated by Formula (1). The aim of Formula (1) is to find all the strips of *O<sup>k</sup>* that are related to *Ok*<sup>+</sup>1; this is equivalent to finding the IDs of the far-left strip and the far-right strip of *Ok*, which are related to *Ok*<sup>+</sup>1.

$$\begin{cases} \mathbf{x}\_s^k = \operatorname{int}(\frac{\mathbf{w}^{k+1} \times (i-1)}{\mathbf{w}^k}) + 1 \\\\ \mathbf{x}\_c^k = \operatorname{int}(\frac{\mathbf{w}^{k+1} \mathbf{x} i}{\mathbf{w}^k}) + 1 \\\\ n = \mathbf{x}\_c^k - \mathbf{x}\_s^k + 1 \end{cases} \tag{1}$$

With the above parameters, the average state value of the objective operating strip can be calculated by Formula (2), where *V* is the mean of the state values of all the related strips of *Ok*:

$$\overline{V} = \sum\_{j=\mathbf{x}\_s^k}^{\mathbf{x}\_r^k} V\_j^k \doteq n \tag{2}$$

Since 3 is the maximum value of the operating strip state, when *V* = 3, that means that the state values of all the related strips are 3 (*done*). This is the condition required to change the strip state of the later operation (from *unready* to *ready*, or from 0 to 1), and the state of *Sk*+<sup>1</sup> *<sup>i</sup>* can then be updated from *unready* to *ready*.

#### *2.3. Sequential Operations Connection*

#### 2.3.1. Field Partitioning

For single-operation tasks, the field can be divided into partitions (sub-fields) according to the number of machinery units, and no optimization is needed within a partition. For multi-operation tasks, a partitioning strategy is also used for the flow-shop working mode. For instance, for the first operation task, all of the machines should work together from one side of the field to the other side, and each machine should occupy several operating strips as a temporary partition (Figure 4). When the machines finish a temporary partition, they move to the next temporary partition. Therefore, in each temporary partition, there will be no collisions between the multiple machinery units.

**Figure 4.** Dynamic partitioning for several sets of machinery units.

In practical work, the partitions for each operation are related to the number of machines, which is different in each operation. Different operations also have different work widths; thus, the partitioning for each operation is dynamic. Figure 4 shows an example of a fieldwork pattern with three sets of agriculture machinery. The entire field area is divided into six static partitions, and each partition is divided into three dynamic sub-partitions (excluding the last partition on the right). Each sub-partition is allocated to a single machine unit.

This study adopted the "first turn skip pattern" (FSP) [28] as the fieldwork pattern, which provides the highest savings in terms of in-field traversing distance among the most common pre-defined fieldwork patterns [29].

#### 2.3.2. Operations Scheduling

The next step of the method includes the calculation of the optimal waiting time between two adjacent operations with a final goal to reduce the overall operating time. As shown in Figure 5, *T*<sup>5</sup> is the overall operating time of the fleet, and for its minimization, it is necessary to shorten the waiting times (*T<sup>k</sup> wait*) between any two adjacent operations.

**Figure 5.** Waiting time and working time of each operation.

To calculate the waiting time, a new parameter, *ck*, was defined to represent the operating capacity per unit time of the *k*th operation:

$$x^k = n^k \cdot v^k \cdot w^k \tag{3}$$

where *nk*, *vk*, and *w<sup>k</sup>* are the number, velocity, and width, respectively, of the machinery unit of the *k*th operation.

Figure 6 shows the optimal waiting time between two adjacent operations. In Figure 6, each color represents an operation, and the horizontal line segment represents the duration within a single static partition. The left oblique line and the right oblique line go through the starting points and ending points of the horizontal line segments. To avoid collisions during operation, the left oblique line of the current operation should not intersect the right oblique line of the previous operation. Therefore, *Tk wait* is determined by the relative relationship of the two slopes of the above two lines and can be calculated using Formula (4):

$$T\_{\text{unift}}^k = \begin{cases} T\_B^{k-1}, & c^{k-1} > c^k\\ \binom{A}{c^{k-1}} - \left(\frac{A}{c^k}\right) - T\_B^{k-1}, & c^{k-1} < c^k \end{cases} \tag{4}$$

where *T<sup>i</sup> <sup>B</sup>* is the duration of a single static partition.

In fact, *ck* is the slope of the left or right oblique line of the *k*th operation.

Due to each operation needing to cover the same area of the field, the total duration (*T<sup>k</sup> work*) of the *k*th operation is calculated as

$$T\_{work}^{k} = A / c^k \tag{5}$$

where *A* is the area of the field.

Finally, the total operating duration of the *k*th operation is then calculated by

**Figure 6.** The relationship between the waiting time and the operating capacity per unit time.

#### *2.4. Collaborative Operating System*

As shown in Figure 7, the cloud-based coordination system has functionalities to manage the field boundaries, operation tasks, agricultural implements, guidance lines, etc. In this study, only A-B guidance lines were used.

(**a**)

**Figure 7.** *Cont*.

**Figure 7.** Part of the software interface of the auto-steering-based collaborative operating system for fleet management (FMCOS): (**a**) Software interface for farmland management and (**b**) Software interface for three-operation coordination.

The more important function of the collaborative operating system is to realize a flow-shop working mode with multiple operations and multiple machinery units. Figure 8 shows the module structure of the collaborative operating system, in which the collaborative module for operations can coordinate all the agricultural machinery within the same operation, and the collaborative module for the fleet can coordinate the agricultural machinery between two sequential operations. Obviously, the fleet has (*N* − 1)-many pairs of combinations of adjacent operations.

**Figure 8.** Module structure of the collaborative operating system for fleet management.

#### *2.5. Auto-Steering System*

As shown in Figure 9, a hydraulic auto-steering system (HuaCe NX200) was used for testing this system. This auto system has the functionality of sending GNSS positions to and downloading guidance lines from the cloud server. Dual-frequency differential signals with an accuracy of ±1.0 cm from the RTK-GNSS base station located at the experimental field were used.

However, NX200 is a commercial product and cannot be modified. For instance, the guidance line (A-B line) cannot be replaced by the operating strip, and the states also cannot be differentiated with colors. Therefore, a T715C tablet (see Figure 9) with an embedded GPS receiver of a single frequency was used, and an Android-based application was developed for better display. The application can download the position data of the NX200 on the same tractor from the server to replace the T715C's own position data to solve the problem of low accuracy. Although the data transform caused some 2 s

delay, the application was able to improve the positioning accuracy of the T715C to ±1.0 cm, which means that the tablet had the same accuracy as the NX200 auto-steering system.

**Figure 9.** The NX200 auto-steering system for a John Deere 1204 tractor and the Android-based application for a T715C tablet.

#### **3. Experiments**

#### *3.1. Simulation Experiments*

A simulation model (Matlab 2018a (MathWorks®)) was developed to validate the state-updating algorithm. The modules of the simulation program are shown in Figure 10. In the model, a strip was simulated as the "pushbutton" of UIControl (see Figure 11), while the operation IDs, the strip IDs, and the values of states were stored in the "UserData" field of the control as a three-dimensional array.

**Figure 10.** Modules of the simulation model of the collaborative operation.

**Figure 11.** Simulation of the state-updating algorithm using Matlab (MathWorks®).

As shown in Figure 12, in the same operation, the MouseClick command was used to change the state from *ready* to *doing* (simulating the machinery operating on the strip), or from *doing* to *done* (simulating the machinery finishing and leaving the operating strip). When the state of *done* appears, the field operation subsequent to the current field operation will check its *unready* strips to determine the state changes by the algorithm.

**Figure 12.** Deterministic finite state automata for the simulation of multi-operation collaboration.

#### *3.2. Field Experiments*

Three experimental setups (see Table 2) were designed to test the FMCOS. These setups, ES1, ES2, and ES3, were designed to test the collaboration functions and performance with one operation, two operations, and three operations, respectively. Three sets of NX200 were mounted on three JD1204 tractors (denoted A, B, and C). The superscript of the tractors' symbols (e.g., A1) represents the sequence of operations.



\* In ES1, tractors A and B worked in the first operation. In ES2, tractors A and B worked in the first operation, and then tractor C worked in the second operation. In ES3, operations 1, 2, and 3 were worked by tractors A, B, and C, respectively.

Figure 13 shows the RTK-GNSS base station, which was the experimental field (located at 40.237789◦ N, 116.581026◦ E) where all experiments took place, and Figure 14 shows the tractors used for tests.

**Figure 13.** Real-Time Kinematic–Global Navigation Satellite System (RTK-GNSS) base station for field experiments and the experimental field.

**Figure 14.** Field experiments scene of (**a**) ES1, (**b**) ES2, and (**c**) ES3.

#### *3.3. Results and Discussion*

#### 3.3.1. Simulation Experiments

Figure 15a presents the scenario of two sets of agricultural machinery for the first operation (*O*1)—which was worked from the middle of the field to both sides—and one set of agricultural machinery for the second operation (*O*2). When clicking on the strip to change the state from *doing* to *ready*, such as *S*<sup>1</sup> <sup>6</sup> <sup>∼</sup> *<sup>S</sup>*<sup>1</sup> <sup>11</sup>, the state of strips related to the second operation, such as *<sup>S</sup>*<sup>2</sup> <sup>6</sup>, changes from *unready* to *ready*, which indicates the appropriate strip for the agricultural machinery of the second operation to work.

**Figure 15.** Simulation result in the flow-shop (**a**) Two-operation and (**b**) Three-operation working modes.

Figure 15b simulates the scene of three sets of machinery for the first operation working simultaneously. Two of them worked from the left or right side to the center of the field, while the other worked from the center to the left side of the field. Similarly, the two sets of machinery for the second operation were able to easily find the *ready* strip for working. The machinery for the third operation can choose the appropriate *ready* strips (such as *S*<sup>3</sup> <sup>11</sup> and *<sup>S</sup>*<sup>3</sup> <sup>12</sup>), and so on.

The results of the above simulations show that the state-updating algorithm can not only well coordinate the state change for multiple machinery within the same operation, but it can also well coordinate two or more operations executed by a fleet of machines.

#### 3.3.2. Field Experiments

Figure 16 depicts screen captures taken during the field experiments of the monitoring interfaces of the collaborative operating system for ES1, ES2, and ES3. ES1 demonstrated the collaborative action of two sets of agricultural machinery (A1 and B1) within the same machinery operation (*O*1). Figure 16a shows the collaborative interface and Figure 17a shows the trajectory of ES1. From the above referenced pictures, we can see that the tractors could get the same base line from the server and generate the same guidance lines in the NX200 terminal. When the tractor entered or left a strip, the strips of the server and the guidance lines of the terminal changed the state. Therefore, operators were able to choose the *ready* strip to work. ES2 demonstrated the collaborative action of three sets of agricultural machinery (A1, B1, and A2) within two machinery operations (*O*<sup>1</sup> and *O*2). Figure 16b shows the collaborative interface, and Figure 17b shows the trajectory of ES2. The experiment showed that the strip state of the second machinery operation could be updated with the operating progress of the first machinery operation. Therefore, the operators within the same operation are able to choose the *ready* strip for work, and the operators of different operations are also notified in order to select the appropriate strips to enter.

**Figure 16.** Monitoring interface of (**a**) ES1 and (**b**) ES2.

**Figure 17.** Trajectory of (**a**) ES1 and (**b**) ES2.

In Figure 17b, in order to show the difference with the second machinery operation, the trajectories of the two tractors (A1 and B1) for the same machinery operation were fused.

As shown in Figure 18, ES3 demonstrated the collaborative action of three tractor operations. The experiment process verified that FMCOS could support the coordination of multiple operations. In fact, three-operation collaboration consists of two two-operation collaborations. Therefore, the performance of ES2 has already proved the feasibility and usability of FMCOS.

**Figure 18.** Monitoring interface of ES3.

The results of the three experiments are shown in Table 3. In ES1, tractors A and B worked in one operation. ES2 consisted of two operations, for which tractors A and B worked in the first operation, then tractor C worked in the second operation. ES3 consisted of three operations, worked by tractors A, B, and C, respectively. Table 3 provides the range of the efficiency improvement for different operating tasks, which have different operation combinations. The maximum efficiency improvement occurs in the experiment of ES3, while the minimum efficiency improvement occurs in the experiment of ES2. In general, the range of efficiency improvement was between 23.05% and 32.72%. It can be concluded that the FMCOS can significantly shorten the waiting time between adjacent operations so as to improve the whole work efficiency.


**Table 3.** Comparison between Non-FMCOS and FMCOS data in regard to efficiency improvement for ES1, ES2, and ES3.

<sup>1</sup> The sum of operating time; <sup>2</sup> Time from the start to the end of the experiment; <sup>3</sup> [(Non−FMCOS)−(FMCOS)] Non−FMCOS <sup>×</sup> 100%.

#### **4. Conclusions**

A cloud-based collaborative operating system for fleet management (FMCOS) was developed by utilizing a flow-shop working mode and by implementing auto-steering systems. The core technology of the system lies in the state-updating algorithm for the operating strips between adjacent machinery operations and the strategy for adjacent operation connection. The system defined the strip states with four classes: *unready*, *ready*, *doing*, and *done*. The FMCOS was developed after testing the algorithm in a simulation model developed in Matlab (MathWorks®). Three field experiments were carried out, namely, one-operation, two-operation, and three-operation experiments, to verify the functionality and performance of FMCOS. The assessment of the performance of FMCOS showed efficiency improvement up to 25.03% and 32.72% when compared with non-FMCOS. The results of the simulation and actual

experiments showed that FMCOS can support in-field fleet management, which can improve the performance of machines in a fleet operating in the same field.

**Author Contributions:** Conceptualization, C.W. and Z.C.; methodology, C.W. and Z.C.; software, Z.C. and L.Y.; validation, Z.C., D.W., B.S., and Y.L.; formal analysis, Z.C.; investigation, C.W. and Z.C.; resources, C.W.; data curation, Z.C., D.W., and B.S.; writing—original draft preparation, C.W. and Z.C.; visualization, Z.C.; supervision, D.D.B.; project administration, C.W.; funding acquisition, C.W., and Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded in part by National Key Research and Development Program of China, grant number 2016YFB0501805. The research was also partly supported by Chinese Universities Scientific Fund, grant number 2018XD003.

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

#### **References**


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

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Energies* Editorial Office E-mail: energies@mdpi.com www.mdpi.com/journal/energies

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18