**1. Introduction**

The optimal allocation and power control of distributed generations (DGs) in electrical power systems has been the focus of researchers in recent years, especially in the context of the smart grid concept. In the project planning phase of DGs installation in the power network, the optimal allocation should determine the system nodes/buses and sizes/powers of the DGs. After the DGs are installed, in the operational phase of the projects, the optimal control of the controllable DGs outputs is of interest. There are different approaches to solve such optimization problem, considering these two parts simultaneously or separately. Besides this question, there is a research challenge regarding the input data resolution used in the optimization process. Generally speaking, the simultaneous approach and higher input data resolution require more computational effort to solve the problem. As far as modeling of the system in the optimization problem is concerned, there are two approaches used in the research. One is based on the usage of an analytical model (a system of equations) of the power network, and the other uses the simulation tool to calculate the objective and constraint functions of the optimization problem. The analytical approach usually means more approximations, and neglecting will be included in the model of the system, decreasing the realistic representation of the system. On the other hand, nowadays, the

**Citation:** Barukˇci´c, M.; Varga, T.; Jerkovi´c Štil, V.; Benši´c, T. Co-Simulation Framework for Optimal Allocation and Power Management of DGs in Power Distribution Networks Based on a Computational Intelligence Techniques. *Electronics* **2021**, *10*, 1648. https://doi.org/10.3390/ electronics10141648

Academic Editor: Cheng Siong Chin

Received: 24 May 2021 Accepted: 8 July 2021 Published: 10 July 2021

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

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

Faculty of Electrical Engineering, Computer Science and Information Technology Osijek, Josip Juraj Strossmayer University of Osijek, 31000 Osijek, Croatia<EMAIL_ADDRESS>(T.V<EMAIL_ADDRESS>(V.J.Š<EMAIL_ADDRESS>(T.B.)

simulation tool for power system analysis ensures less neglecting, resulting in more realistic modeling of the system. The results obtained using more realistic model in the optimization process are more reliable for practical implementation. These two research issues in this topic are the research subject of this paper. Besides the single or multi-objective approach, the different objective functions as well as the problem constraints of the optimization problem are used in research studies of optimal allocation and the control of DGs. A brief overview of the literature on the topic concerning the above-mentioned issues and research challenges is presented below in the section.

In [1] the authors used the fuzzy system to aggregate the multi-objective problem into single-objective optimization, considering the economic and environmental objective functions and technical constraints. The objective functions were calculated from analytical expressions, and node sensitivity analysis (gradient-based optimization) was used to find the optimal DG allocations. The optimization assumed constant (static) load powers in the power network. The load powers dependent on a bus voltage were considered in [2] to find the optimal allocation of DGs. However, the nominal loads were taken to be constant—not changing in time. The single-objective objective function value and technical constraints were used and solved by the metaheuristic optimization method (Harmony Search Algorithm—HAS). The calculation of the objective function values was coded in the programming environment in which the optimization method was implemented. The multiobjective optimization problem of optimal allocation and control of Battery Energy Storage Systems (BESS) is presented in [3]. Particle Swarm Optimization (PSO), which belongs to the metaheuristic techniques, was used to solve the problem with two objectives: power losses and total power of installed BESS. Optimal allocation was performed with constant nominal load values. After the optimal BESS allocations were found, the optimal charging/discharging control of BESSs with the hourly resolution was determined, using the proposed analytical method. In [4], the single optimization problem considering the DG penetration level as objective and total harmonic distortion (THD) as the problem constraint were investigated again with static load values. The problem was solved by metaheuristic optimization techniques, namely Genetic Algorithm (GA) and Differential Evolution (DE). In addition, in [3,4], objective function values were coded together with the optimization method. In [5], the multiobjective optimization of the optimal power generation of Virtual Power Plant (VPP, including DG and EV (Electrical Vehicle) charging stations) power generation was solved, using the PSO method. The daily load shape with the hourly resolution was used in the simulations. Two-objective optimization, considering the operational cost and which are the pollutants emission, was presented in [6]. The problem was solved using the Ant Lion Optimizer (ALO) metaheuristic method, and variable loads on a daily level with hourly resolution were used. The optimization procedure, as well as the problem, were coded together inside the same programming tool. The single objective optimization dealing with losses minimization and DG penetration level maximization was presented in [7] to find optimal allocation of DGs. The daily profile of load, Photovoltaic (PV), and Wind (W) generations DG units with minute resolutions for 12 typical days were used in the research. The objective function values were calculated by coding the problem in a programming tool, and it is not clear what optimization method was used here. The single objective optimization problem of optimal DG allocation was solved by applying the metaheuristic optimizer, Grey Wolf Optimizer (GWO) in [8]. The optimization was performed for a constant load value, and methods for power flow calculations and optimization itself were coded with the same computational tool. In [9], single-objective optimization, considering the optimal allocation of DG for reactive power control, was solved by using GA. The co-simulation approach was used here, employing the external power system simulator to calculate the objective function value. As in most literature considering the topic, the optimization was performed for constant load value. The constant loads were considered in [10] during single-objective optimization, aiming to find the optimal DG allocations with power losses minimization. The objective function values were calculated based on the backward/forward sweep power flow method which

was coded alongside the optimization method in the same programming environment. Different metaheuristic techniques, HAS, Artificial Bees Colony (ABC), and PSO, were used to solve the optimization problem. The research [11] presented a co-simulation approach to solve the optimal DG allocation problem considering one objective function and constant loads. The single objective optimization of DG allocation was solved by using the Sensitivity Analysis (SA) method in [12]. The optimization method and calculation of the objective function value were implemented in the same programming tool. The yearly profile (at hourly data resolution) of wind DG production was considered, but it is not clear if the load profile was used too or if the loads were assumed to be fixed. In [13] the multiobjective optimal DG allocation problem was solved by using the metaheuristic optimization method Ant Lion Optimizer (ALO). The optimal allocation of DGs, BESS, and reactive power control devices was found, considering constant loads. The paper [14] dealt with optimal DG and BESS allocations, solving the problem by using mixed-integer conic programming (MICP) as an optimization technique. The optimal allocation problem was defined in a form of single-objective optimization. The optimization problem was modeled in the specific modeling tool, and an existing external optimization tool was interfaced to the model to find the optimal solution. The optimization problem considered load and DG production profiles with hourly resolution at the yearly level. However, the clustering technique was used to generate typical 48 profile patterns to decrease the dimensionality of the problem. The authors in [15] used PSO to solve the single optimization problem of optimal DG allocation. The optimization process considered constant load and power flow calculation, as well as the optimization algorithm, and was coded in the same programming tool. In [16], the hybridization of two metaheuristic methods, PSO-SFL (Shuffled Frog Leap), was used to solve the single optimization problem of optimal DG allocation. The constant load values were considered and power flow calculations were implemented in the programming tool used for the optimization method performing. The research presented in [17] solved the optimal allocation of DGs, Shunt Capacitors (SC), and Electrical Vehicle (EV) charging stations by using the Grasshopper Optimization Algorithm (GOA), which belongs to the class of metaheuristic methods. The proposed procedure solved the singleobjective optimization problem (with four objectives aggregated into one objective function) in two separated steps considering the optimization of DG and SC allocations separately for the optimization of the EV charging station allocation. The constant load values were used during the solving of the optimization problem, and the impact of the load and DG production changes was investigated once the optimal DG and SC allocations were determined. The power flow analysis was performed in the same programming tool, which was used for the optimization method implementation. A hybrid metaheuristic method, GA-PSO, was applied in [18] to solve a single-objective problem with the aggregated objective function consisting of three parts. The load values were assumed to be constant, and the power flow calculation was implemented in the optimization method objective function value calculation. The authors in [19] applied GA optimization to solve the single optimization problem of optimal DG and BESS allocations, considering the daily load shapes with hourly resolution. The objective function calculation was coded inside the optimization procedure in the PYTHON programming environment. In [20], the single-objective optimization problem of optimal DG allocation was solved by the DE optimization algorithm. The optimal DG allocation, as well as the DG power factor, were decision variables of the optimization problem. The load and DG production profiles were considered during the optimization. The daily profiles with hourly resolution were used, and these profiles were obtained by averaging values from the seasonal profiles. The DIgSILENT simulation tool was employed for the power flow calculation based on which the objective function value was calculated. Although not clearly stated, if the optimization method was implemented inside (built-in) the DIgSILENT or in some external programming tool, it seems that the co-simulation approach was used here.

Based on the above-given brief overview of research studies dealing with optimal DG allocation, the description of the reviewed literature can be summarized as follows:

The above literature review yields three main issues considering the optimal allocation of DGs:

These previously mentioned are detected as open research questions, challenges, and gaps in the topic of optimal DG allocation. The presented research aims to decrease the research gaps and make a contribution to the topic through the following:

The existing literature includes the above-listed aspects (some or all). The optimal DG allocation problem considered in this paper can be summarized as follows:

The approach of simultaneously considering the location, size, load, and DG production profiles during the optimization used in this research study is similar and with similar aims as that presented in [20], which was one of the inspirations for this research.

The rest of the paper is organized as follow: the optimization problem formulation is given in Section 2; Section 3 describes the applied research methodology and proposed framework, including a brief overview of the used specific simulation tools; the results of different scenarios aiming to validate the proposed framework applied on a test power distribution network are presented in Section 4; the discussion about the obtained results concerning the stated research question regarding the input data resolution impact is given in Section 5; and at the end, some general remarks/conclusions are presented in Section 6.

#### **2. The Optimization Problem of the DGs Allocation**

The optimization problem is defined to address the research gaps mentioned above, in the previous section. In the literature, different objectives are considered; here, two objective functions that are important from the point of view of the power distribution system operator and the owners of the DGs are used. As one of the main interests of the power distribution system operator is decreasing the losses, the first objective function used in optimization is the active energy losses *Wloss*. Generating as much energy from DGs as possible is the main interest of the DG owners to shorten the investment payback period and increase the profit. This leads to increasing the penetration level of DGs in the power distribution system. Such an objective is formulated here in the form of the total exchange of the total apparent energy *WS* in the power coupling point of the distribution network on the upstream network. Because these two objective functions are conflicted, the multiobjective (two-objective) optimization approach is used in the research. The technical constraints regarding the nodal voltage range and line current limits as well as the box constraints of the decision variable ranges are applied in the optimization problem.

The mathematical notation of the multiobjective optimization problem, including the previously described objective functions and constraints, is as follows:

$$\begin{array}{l} F = \left[ W\_{loss}(\mathbf{x}), W\_S(\mathbf{x}) \right] \to \min \\ \text{subject to box constraints:} \mathbf{x} \in \{ \mathbf{x}\_{lb}, \mathbf{x}\_{ub} \} \\ \text{subject to inequality constraints:} \\ V\_{min} \le V\_{i\varepsilon} \le V\_{max}, I\_{k\varepsilon} \le I\_{k,max} \le \varepsilon \\ \text{with decision variables vector: } \mathbf{x}\_{\varepsilon} \end{array} \tag{1}$$

with the following notations: *F*—two objective problem function consisting of *Wloss*, yearly energy losses in the network, and *WS*, yearly exchange of the apparent energy between a network with DG and upstream system; *xlb*—lower bounds of the decision variable values; *xub*—upper bounds of the decision variable values; *Vmin*—lower bounds of the nodal voltage value; *Vmax*—upper bounds of the nodal voltage value; *Vi*,*e*—calculated nodal voltage in the *i*-th network node; *Ik*,*max*—maximum allowed currents in the *k*-th network line; and *Ik*,*e*—calculated current in the *k*-th network line.

The objective functions represent energies over the timespan and for *N* time segments *ti* are calculated from active power losses (*Ploss*), active (*Pexc*) and reactive (*Qexc*) powers exchanged with the upstream network as follows:

$$\mathcal{W}\_{\text{loss}} = \frac{\mathcal{W}\_{\text{loss},i} \cdot t\_i}{\sqrt{\left(\sum\_{i=1}^{N} P\_{\text{exc},i} \cdot t\_i\right)^2 + \left(\sum\_{i=1}^{N} Q\_{\text{exc},i} \cdot t\_i\right)^2}} \tag{2}$$

Optimization Problem (1) has a two-objective function that consists of two objectives: yearly energy losses in the network (*Wloss*) and yearly exchange of the apparent energy between a network with DG and upstream system (*WS*). The problem constraints in (1) are related to the decision variable ranges (box constraints) and the network operational constraints. The box constraints represent ranges of the decision variable values. The operational constraints are related to the standardized nodal voltage ranges (*Vmin* − *Vmax*; usually nodal voltage limits are given in range ±5% or ±10% depending on the relevant standard) and rated currents of the network lines (*Ik*,*max*).

Optimization Problem (1) is solved by using the Pareto dominance definitions resulting from the solution set known as the Pareto set [21]. Except for the multiobjective optimization, the two single objective optimizations considering each of the objective functions individually are performed also to check if the multiobjective approach can find reliable edges of the Pareto set. The three different problem setups are used in the research, each of them resulting in a different number of decision variables. This being the case, the contents of the decision variables vector is detailed later in the text in Section 3, describing the proposed procedures. The co-simulation setup of the distribution power system simulation and tools of the computational intelligence methods are used to solve the black-box optimization problem model.

The simulation tool for electrical power system simulation is used to calculate the objective function values. Both DGs types with controllable and non-controllable primary energy sources are considered in the distribution system model. The non-controllable DGs used in the research are Photovoltaic (PV) and wind (WD) plants; as controllable DGs, Biogas (BG) plants are used in the distribution power network. For non-controllable energy sources, the production profiles (the DGs outputs) are involved in the model; for

the controllable source, the DGs outputs are subject to optimization. In addition, the load profiles (costumers load shapes) are considered in the optimization. The base case in the simulations supposes DG production and load profiles on a yearly level with hourly resolution, i.e., with 8760 data.

#### *2.1. A Brief Overview of the Used Tools in Co-Simulation*

The optimization tool applying the metaheuristic optimization technique MIDACO solver (Mixed Integer Distributed Ant Colony Optimization) [22] is used to solve the optimization problem. The advantage of this tool is its applicability in a general case of complex optimization problems, including continuous (linear (LP) and nonlinear (NLP)), integer (discrete) (IP), and mixed-integer (MINLP) problems. This tool is based on the Ant Colony Optimization (ACO) [23] which belongs to the metaheuristic methods. The MIDACO can handle single as well as multiobjective optimization.

The artificial neural network (ANN) is applied in the optimization problem, aiming to decrease the number of the decision variables, reducing the dimensionality of the problem. The Multilayer Perceptron (MPL) ANN with one hidden layer is used here. The ANN is modeled by usage of the TensorFlow tool [24] by applying the Keras API [25].

The OpenDSS simulation tool [26] is used in the study with the purpose of solving the power flow calculations of the model of the distribution network and obtain the objective function values needed for the performing optimization procedure. The usage of the power network simulation tool enables including more details in the network model, decreasing the approximations and neglecting in the model, compared to an analytically defined objective function. This ensures more realistic power distribution network modeling, resulting in more reliable results obtained through the optimization.

As mentioned before, all three of these computational tools are implemented in the Python programming environment and employed in the co-simulation setup to solve the black-box optimization problem.

#### **3. Co-simulation Framework for DG Allocation and Power Management Optimization**

The research methodology used in the study is based on the implementation of different scenarios for optimization and the investigation of the input data resolution impact on the solution of the optimization problem. The optimization scenarios are related to variable or fixed power factors of each DG and directly determine the number of decision variables in the optimization problem. The three cases described in Table 1 are used as optimization scenarios.

**Table 1.** The different optimization scenarios.

As can be seen in Table 1, the values of the DG power factor can be directly included as a decision variable, and in that case, the fixed value of the power factor is optimized and does not change with DG power output changes. The other scenario investigates application of the ANN to estimate the optimal DG power factor and in this case, the DG power factor changes over time. For all scenarios, the output power of controllable DG is determined by the ANN. These three basic cases will result in different types and numbers of decision variables in the optimization problem. In Table 2, an overview of the decision variables that occur in the described optimization scenarios is presented.

**Table 2.** The optimization decision variables for different scenarios.

In all scenarios, the consumers' load shapes and energy source profiles are used to model the time variability of the primary energy sources and consumption.

In Figures 1–3, the proposed frameworks according to the previously described scenarios are shown. Application of the ANN to estimate the optimal DG power output and DG power factor values is one of the main contributions of this research. The purpose of the ANN is to decrease the number of decision variables in the optimization problem. Without the proposed usage of the ANN for each controllable (in some scenarios, also for uncontrollable (Table 2)) source, the number of decision variables will be equal to the input data resolution, e.g., for yearly input data with the hourly resolution it will be an additional 8760 decision variables per DG for scenarios Opt 1, i.e., twice for scenario Opt 2 or four times for scenario Opt 3 (Table 2).

Based on the scenarios overview given in Table 1 and Figures 1–3, the similarity and differences between the optimization models Opt 1–Opt 3 can be summarized as follows. For all optimization scenarios, uncontrollable DG outputs are considered variable in time and defined by energy source profiles. The controllable DG output is managed by the ANN in all scenarios. The differences among the scenarios are related to the power factor variability for controllable and uncontrollable DGs. In Opt 1, the fixed DGs power factors of both controllable and uncontrollable DGs are optimized. The fixed and ANN managed power factors of the controllable and uncontrollable DGs, respectively, are optimized in scenario Opt 2. In Opt 3, modeling the power factors of all DGs (controllable and uncontrollable) is managed by ANN.

The details about the ANN inputs and outputs are given in Figure 4, and the ANN parameters (weights and biases) are optimized by the MIDACO solver simultaneously with DGs allocation optimization. The simple Multilayer Perceptron (MLP) ANN consists of one hidden and output layer. Because the usage of the ANN has the purpose of significantly decreasing the number of decision variables, this simple ANN configuration is implemented in the framework. The purpose of ANN is to significantly decrease the number of problem decision variables. If the ANN is not used, the number of decision variables would be increased for the number of the input data (depending on the data resolution, e.g., in case of yearly data with hourly resolution the optimization problem would have 8760 data only for the DG power management). The ANN takes as inputs the *i*-th data from each load shape and each uncontrollable DG output profile. Training of the ANN was performed in this way (without application of standard built-in ANN optimizer presents in the used ANN tool) because the training process, in this case, is slightly different than ordinary. Usually, when ANN is used for regression purposes, the difference between the target and ANN estimated values is objective in the ANN optimization process (training). However, this is not the case here; the ANN target outputs are not known in advance, as they need to be determined during the optimization of the whole problem defined by (1). Such configuration of the proposed framework and purpose of the ANN prevents performing common ANN training procedures.

**Figure 1.** The framework used for Opt 1 optimization scenarios.

**Figure 2.** The framework used for Opt 2 optimization scenarios.

**Figure 3.** The framework used for Opt 3 optimization scenarios.

**Figure 4.** (**a**) Overview of data types used by the ANN; (**b**) schematic structure of the ANN.

The two different workflows are applied to investigate the possible improvement of the solution quality for scenarios Opt 2 and Opt 3 (Figure 5). The difference between the workflows is about the initial solution that the optimization algorithm starts with. In both workflows (WF1 and WF2), the initial solution is randomly generated in Opt 1 optimization scenarios. In workflow W1 the initial solution is also randomly chosen for scenarios Opt 2 and Opt 3, while in workflow WF2 the starting optimization point for Opt 2 and Opt 3 is the solution of solved Opt 1. In the case of WF2 workflow, the locations, and sizes of all DGs and trained ANN for controlling the power output of controllable DG are the solution from Opt 1, and only parameters of the ANN are used to optimize the DG power factors. The purpose of the WF2 workflow is to check if the optimal solution obtained by solving Opt 1 can be additionally improved by managing the DGs power factor.

The above-described optimizations and workflows are performed with input data on yearly level with hourly resolution, i.e., with 8760 input data modeling changes of primary energy source intensity and loads over a year. After the research done according to the workflows, the final workflow and optimization model will be proposed based on a comparison of the obtained solutions for different optimization scenarios. In the continuation of the proposed procedure, the previously chosen optimization scenario and workflow will be repeated with the different decreased resolutions of input data to research the impact of data resolution on the optimization problem solution. The two decreased input data resolutions with 12 and 365 input data are used here. These input data are obtained by averaging the basic input data (hourly resolution, 8760 input data) on monthly (12 input data) and daily (365 input data) levels. The schematic overview of this part of the research methodology is given in Figure 6.

**Figure 5.** The workflows used in the research study, (**a**) workflow uses randomly generated input data for single and multi-objective optimizations; (**b**) workflow uses solution of the multi-objective optimization as input to single-objective optimizations.

**Figure 6.** Application of input data with different resolutions in the proposed procedure.

At the end of the research methodology used in this study, the chosen optimization procedure optimized with hourly input data was applied, using input data with a resolution of 15 min (four times higher resolution than one used in optimization procedure) steps, i.e., with 35,040 input data. The load shapes are one of the inputs into the ANN (Figure 4). Usually, the load shapes are forecasted, based on the historical consumption data, and consequently include more or less uncertainty. Due to this uncertainty, the later application of the load shape can give a load value that is different from the real load value. This procedure part is performed with two scenarios: one supposing an unknown real load value at the moment and the other with a known specific load value at the specific time step. If the real load value is not known/measured, the load value is estimated according to the load shapes used in the optimization process. The previously described is visualized in Figure 7. The purpose of this is to investigate the robustness of the proposed framework. In addition, with the development of the smart grid concept, it is expected more and more the usage of smart energy meters at the costumers' point of common coupling. These conditions allow obtaining the load values with more accuracy as input data during the implementation of the proposed framework.

**Figure 7.** Scenarios for the procedure application for 15 min resolution of input data.

#### **4. Application of the Framework on Test Distribution Power Systems—Case Study**

The very well-known and often used IEEE 13 node IEEE 37 node test feeders [27] are used as a case study to demonstrate the application of the proposed framework for optimal DG allocation and power management. These distribution power systems are chosen because they represent the most general examples of the power networks. Some features of the used power networks are an unbalanced system (one, two, and three-phase lines and loads), different load models (constant power, constant impedance, constant current ...), voltage regulators, load in wye (star), and delta connection configurations and two voltage levels (4.16 kV and 0.48 kV). All details about the used test systems can be found in [27]. The next modification in the original IEEE 13 node and IEEE 37 node bus test feeders is made to adapt the systems for performed research: the capacitor banks are turned off, the taps of the voltage regulators are reset to the middle position before the simulation for each

possible solution in the optimization is started, and three different load shapes are assigned to the consumers according to Table 3. The proposed optimization procedure is performed on a desktop PC with Intel i7-10700 CPU 2.90 GHz, 8 Cores, RAM 16 GB. The versions of the used software are PYTHON 3.8, MIDACO 5.0, OpenDSS <IP_ADDRESS>, and Tensorflow 2.3.1. Data about the optimization parameters used in the optimizations are as follows: number of function evaluation in MIDACO is 30,000, number of ants and kernels in MIDACO are 250 and 15, respectively. The computational times are 1.1 and 6 s per iteration for IEEE 13 and 37 bus test networks, respectively. The total computational times for these two tested networks are about 9 and 22 h for IEEE 13 and 37 bus networks, respectively.

**Table 3.** The load profile of costumers used in the test networks.

#### *4.1. Input Data Preparation*

The input data regarding the load shapes and PV and wind plant power profiles are obtained by usage of the existing tools for prediction/forecasting load consumption and PV and wind DG production. The computer tool Load Profile Generator (LPG) (https://www.loadprofilegenerator.de/ (accessed on 8 April 2021)) [28] is used to synthesize the three different load profiles used in the simulations. The built-in (LPG) load shape types, namely "H01 in HT 14", "H01 in HT 11", and "H01 in HT 07", are used to generate load profiles LP 01, LP 02, and LP 033, respectively. The normalized load shapes are given in Figure 8a–c by showing example daily profiles in Figure 8e–f.

The online platform (tool) named "Renewables.ninja." (Available online: https:// www.renewables.ninja/ (accessed on 8 April 2021) based on research presented in [29,30] is used to produce the output profile of PV and wind DGs. The generated normalized DG production profiles are shown in Figure 9.

The above-described input data, the yearly load, and DG production profiles with the hourly resolution are generated as basic input data cases. The input data with monthly and daily resolutions used as input into part of conducted research study shown in Figure 6 are generated from the base input data by averaging data on monthly and daily levels, respectively. These input data with the decreased resolution are presented in Figures 10–13.

The input data used at the end of the conducted research shown in Figure 7 are produced by adding randomly generated noise ( *X* according to uniform distribution U) to the base data for each time step (1 h) four times. The repeated random number generation four times for each hour produces load and DG output profiles with a resolution of 15 min (1 h = 4 × 15 min). The noise range is set in the range ±*z*% of the hourly value. The mathematical formulation of the 15 min resolution data is the following:

$$LF\_{15\text{min}} = LF\_{\text{h}} + X \sim \mathcal{U}(-z \cdot LF\_{\text{h}}, z \cdot LF\_{\text{h}}) \tag{3}$$

In this case, it is important to highlight that hourly data are not averages of 15 min data. Due to higher data resolution, the 15 min resolution data are not visualized for the whole year, but only for the first day in a year (in comparison with basic input data (hourly resolution)) as can be seen in Figures 14 and 15 (for example, for ±20% uncertainty used in the 15 min data generation).

**Figure 8.** Input data—load shapes: (**a**) LP 01-year; (**b**) LP 02-year; (**c**) LP 03-year; (**d**) LP 01-day; (**e**) LP 02-day; (**f**) LP 03-day.

**Figure 9.** Input data—load shapes: (**a**) PV plant output; (**b**) wind plant output.

**Figure 10.** Load shapes with daily resolution: (**a**) LP 01; (**b**) LP 02; (**c**) LP 03.

**Figure 11.** Load shapes with monthly resolution: (**a**) LP 01; (**b**) LP 02; (**c**) LP 03.

**Figure 12.** Production profile of DGs with daily resolution: (**a**) PV DG; (**b**) wind DG.

**Figure 13.** Production profile of DGs with monthly resolution: (**a**) PV DG; (**b**) wind DG.

**Figure 14.** Load profiles with 15 min resolution generated (according to (3)) from hourly resolution (for example, for the first day in a year): (**a**) LP 01; (**b**) LP 02; (**c**) LP 03.

**Figure 15.** DG output profiles with 15 min resolution generated (according to (3)) from hourly resolution (for example, for the first day in a year): (**a**) PV DG; (**b**) Wind DG.

In Sections 4.2–4.6, the results of the framework implementation are presented.

#### *4.2. Results for Different Optimization Models—Workflow WF1*

In this subsection, the results obtained according to workflow W1 (Figure 5a) for the optimization models proposed in Figures 1–3 and Table 1 are presented. Figures 16–18 visualize solutions of optimization models Opt 1 (Figure 1), Opt 2 (Figure 2) and Opt 3 (Figure 3), respectively, showing the Pareto front of the solved optimization Problem (1).

In Tables 4–6, the numerical values for the solutions obtained for all three proposed optimization models (Opt 1, Opt 2 and Opt 3) applied on multi-objective and single objective optimizations are given. For multiobjective optimization, the values for the Pareto edges (for solutions giving the lowest energy losses and energy exchange) are shown in these Tables.

**Figure 16.** The Pareto front of the solved multi-objective optimization problem for Opt 1 (Figure 1).

**Figure 17.** The Pareto front of the solved multi-objective optimization problem for Opt 2 (Figure 2).

**Figure 18.** The Pareto front of the solved multi-objective optimization problem for Opt 3 (Figure 3).

**Table 4.** Results for optimal solutions of Opt 1 optimization model for different objective functions (W1 workflow).

**Table 5.** Results for optimal solutions of Opt 2 optimization model for different objective functions (W1 workflow).

**Table 6.** Results for optimal solutions of Opt 3 optimization model for different objective functions (W1 workflow).

*4.3. Results for Different Optimization Models—Workflow WF2*

The results obtained according to workflow WF2 (Figure 5b) for the optimization models Opt 2 and Opt 3 are presented in this subsection. In WF2, the workflow of the optimal allocations of DGs and trained ANN obtained as a solution of Opt 1 model are used as the initial (starting) solution for the optimization of the Opt 2 and Opt 3 models. In this case, the problem solutions are trained ANNs for controlling DGs power factors. In this workflow, only a single-objective optimization approach is applied based on the results obtained in the previous subsection. As can be seen from the results obtained in workflow WF1, the single-objective optimization finds the Pareto front edges of better quality. The results obtained according to workflow WF2 are shown in Tables 7 and 8 for optimization models Opt 2 and Opt 3, respectively.

**Table 7.** Results for optimal solutions of Opt 2 optimization model for different objective functions (W2 workflow).

**Table 8.** Results for optimal solutions of Opt 3 optimization model for different objective functions (W2 workflow).

#### *4.4. Comparison of Optimal DG Allocations for Different Optimization Models and Workflows*

This subsection presents the optimal allocations of DGs obtained according to different optimization models and workflows used in the research study. The DG locations (bus) in the network and nominal DG apparent powers (*Sn*) are presented in Table 9 to give a comparable overview of the obtained solutions of the optimization problem, according to the different research scenarios used in the study.

**Table 9.** Results comparison for optimization models and workflows in case of single-objective optimizations.

#### *4.5. Impact of Input Data Resolution on the Optimization Problem Solution*

The previous procedures used in the presented research are analyzed to propose a suitable optimization model for the next steps in the study. Based on the results presented in the three previous subsections, the optimization model Opt 1 is chosen as the one with the best results (more details about the choice are given in the Discussion section below) for the application in the rest of the study. In this subsection, the results of the performing research step presented in Figure 6 are presented. The purpose of this step is to investigate the impact of the input data resolution on the solution of the optimization problem. The solutions for yearly input data with daily and monthly resolutions are compared against the solutions for the data with hourly resolution (as the base case). The input data with daily and monthly resolutions are generated as described in Section 4.1. Table 10 shows objective function values without installed DGs for different resolutions of input data. In Tables 11 and 12, an overview of solution comparisons presented in Figure 6 are shown for optimization model Opt 1 and single optimization problems (with *Wloss* and *WS* objectives separately).

**Table 10.** The objective function values without installed DGs for the used input data resolutions.

**Table 12.** Impact of the input data resolution (Opt 1 optimization problem with *WS* objective).

#### *4.6. Robustness of the Proposed Optimization Model*

Based on the presented results, the optimization model Opt 1 is proposed as a framework for optimal DG allocation and power management. Let us recall that Opt 1 model results with solutions of the optimal allocation of DG, constant optimal power factors of DGs, and trained ANN for power management of the controllable DGs (Figure 1). After the optimal solution with hourly input data is found, the optimal allocation of uncontrollable DGs (PV and wind plants) and the trained ANN is implemented in the simulation model to drive the output power of the controllable DG (the bio-gas plant) with input data with higher resolution than the one used during the optimization process. The 15 min resolution data are used here, applying two scenarios. The ANN inputs are normalized PV and wind plant power outputs and load shape values. The 15 min PV and wind plant outputs, as well as the load shapes, are generated according to the procedure given in Section 4.1. Both scenarios assume the known real PV and wind plant power outputs. The one scenario (S1) assumes unknown real load shape values, in this case, the load shape ANN inputs are defined, according to the given load shapes (used during the optimization) and different from the real load shape value. In this scenario, for each four real load data in the load shape, the load inputs to ANN are the same (according to the given hourly load shape). The other scenario (S2) assumes known load shape values are used as ANN inputs. Both scenarios S1 and S2 are visualized in Figure 7. The results for these two scenarios for ±20% uncertainty (according to (3)) of load shapes and the DG profile are given in Table 13.

**Table 13.** Results of the proposed model application for 15 min input data for Opt 1 model— uncertainty of 20%.

#### *4.7. Implementation on Middle Sized Distribution Network—IEEE 37 Node Test Feeder*

In this subsection, the proposed procedure is applied on the IEEE 37 Node Test Feeder as an example of the distribution network of middle size. Based on the results and analysis of the procedure implementation on the IEEE 13 Node Test Feeder (given in previous subsections), the optimization model Opt 1, workflow WF 1, yearly profiles with hourly data and single objective ( *Wloss*) optimization are applied on the IEEE 37 network. The obtained results are presented in Tables 14 and 15. The two optimization processes are performed: one for allocation of a total of three DGs (1xPV DG, 1xW DG, 1xBG DG), and the second for a total of six DGs (2xPV DG, 2xW DG, 2xBG DG).

**Table 14.** Optimal allocation of DGs in the IEEE 37 Node Test Feeder.

**Table 15.** Decreasing the yearly power losses by optimal allocation and power management of DGs in the IEEE 37 Node Test Feeder.


**Acknowledgements**

the SQI was evaluated based on the parameters presented here, to which information about the drainage was added (considering the texture, the shape of the terrain, and the average annual precipitation), this did not alter the trend of the results, and some conclusions can be drawn about the consequences of different land uses on the quality of the soil and, in turn, on

The obtained results showed that the higher the OM, the higher the rate of infiltration and the water retention in the soil, apparently with little influence of soil chemical constituents on aquatic ecosystems, since the water that solubilizes or suspends them is mostly retained in the soil. Only in situations where land use practices implies leaving the soil bare during the beginning of the rainy season, it will be possible to verify this type of consequences through

Agricultural land use has shown to have strong negative effects on water quality, stream habitat, and degradation of riparian vegetation, ultimately resulting in fish assemblages' impoverishment and clearly benefiting non-native species, which thrive under altered conditions [95]. The most negative effects were associated with intensive, heavily irrigated, fertilized, and pastured agricultural systems, mostly represented at the basin scale by olive groves, irrigated crops, and pastures. Conversely, agroforestry (montado) results emphasize the potential con-

Since Mediterranean rivers exhibit high levels of fish fauna endemicity, human impacts on these systems have the potential to extirpate native species and reduce local, regional, and global native biodiversity [96]. It should be further highlighted that considering the forecasted climate changes and their possible joint effects with land use changes, far reaching effects are likely to occur on ecological communities in Mediterranean regions in the future [97, 98].

Regarding soil, local agricultural intensity did not prove to be a threat to the integrity and quality of the soil, seeming to ensure the sustainability of the local uses and practices, contrary to what is usually found. The intensification of agricultural systems, by means of a high consumption of water or energy, can be carefully planned, thus preventing soil degradation through the known threats defined by Thematic Strategy for Soil Protection [99], namely, decline in organic matter, compaction, floods, soil erosion, salinization, contamination, landslides, and sealing. The careful planning and execution of agricultural practices that intensify the production systems (but not involving degradation) are the recommendations of several

These are preliminary findings based on a case study, and more detailed research is further required to substantiate the results and assess the direct relationship between soil and the aquatic ecosystem, namely, by considering more soil characteristics and diversity, covering larger spatial and temporal scales, and considering climate data. This would allow a better understanding of the complex pathways underpinning the interaction among the processes

tribution of this stable production system to biodiversity conservation.

the adjacent aquatic ecosystems.

42 Land Use - Assessing the Past, Envisioning the Future

the surface runoff [94].

**4. Conclusions**

authors [88, 100, 101].

and factors involved.

This work is funded by National Funds through FCT - Foundation for Science and Technology, I.P., within the framework of the project UIDMulti044632016. We thank the Institute of Mediterranean Agrarian and Environmental Sciences (ICAAM) and the Center for Interdisciplinary Development and Research on Environment, Applied Management and Space (DREAMS). Special thanks to Maria José Barão for laboratory analyzes of water samples. We are also grateful to all the colleagues and volunteers for their valuable help in the fieldwork.


**Nitrogen and Cobalt Co-Coped Carbon Materials Derived from Biomass Chitin as High-Performance Electrocatalyst for Aluminum-Air Batteries**

#### **Mi Wang 1, Jian Ma 2, Haoqi Yang 2, Guolong Lu 2,\*, Shuchen Yang 1,\* and Zhiyong Chang 2**

Received: 21 October 2019; Accepted: 10 November 2019; Published: 14 November 2019 -

**Abstract:** Development of convenient, economic electrocatalysts for oxygen reduction reaction (ORR) in alkaline medium is of great significance to practical applications of aluminum-air batteries. Herein, a biomass chitin-derived carbon material with high ORR activities has been prepared and applied as electrocatalysts in Al-air batteries. The obtained cobalt, nitrogen co-doped carbon material (CoNC) exhibits the positive onset potential 0.86 V vs. RHE (reversible hydrogen electrode) and high-limiting current density 5.94 mA cm−2. Additionally, the durability of the CoNC material in alkaline electrolyte shows better stability when compared to the commercial Pt/C catalyst. Furthermore, the Al-air battery using CoNC as an air cathode catalyst provides the power density of 32.24 mW cm−2 and remains the constant discharge voltage of 1.17 V at 20 mA cm−2. This work not only provides a facile method to synthesize low-cost and efficient ORR electrocatalysts for Al-air batteries, but also paves a new way to explore and utilize high-valued biomass materials.

**Keywords:** electrocatalyst; oxygen reduction reaction; Al-air battery; biomass; nitrogen-doped carbon

#### **1. Introduction**

Environmental deterioration and the depletion of finite resources stimulate research to seek and develop clean energy storage and conversion devices for the future. Metal (Al, Li, Zn, etc.)-air batteries [1–3] are considered the most promising candidates to replace the current fossil fuels due to their low carbon emissions, high theoretical specific capacity and energy density, etc. Among the batteries listed above, the Al-air battery is most abundant in the Earth's crust and possesses a high energy density (8.10 Wh g−1) [4], which gives the Al-air battery considerable potential to be widely used; however, its poor electrochemical performance constrains its practical application.

Oxygen reduction reaction (ORR) is a key reaction affecting the discharged performance of energy storage devices such as fuel cells and metal-air batteries [5–7]. In order to enhance the discharged performance of these devices, a platinum-based electrocatalyst is commonly applied to overcome the sluggish ORR kinetics and accelerate the reaction process [8–11]. Despite the fact that platinum (Pt)-based catalysts are considered the most efficient ORR electrocatalyst, their expensive price and poor durability has limited their practical application [12–14]. Thus, great efforts have been devoted to exploring low-cost and high-performance carbon-based electrocatalysts. Among various carbon materials, graphene [15], carbon nanotubes (CNTs) [16] and activated carbon [17] are most frequently used in the preparation of carbon-based electrocatalysts due to their excellent physical and chemical properties. Although these carbon materials have shown remarkable catalytic activity, their cost still

limits their application. Therefore, it is important to seek another substitute for synthesizing a low-cost, carbon-based ORR electrocatalyst.

Biomass materials are the ideal carbon precursor for synthesizing ORR catalysts due to their renewability and large-scale and low-cost production [18–25]. Recently, several studies reported biomass-derived ORR catalysts with efficient catalytic activities. Li and co-workers [26] reported a low-cost method for transforming woody biomass into a single-atom, dispersed ORR catalyst. Wu et al. [27] synthesized kelp-derived, multilayer-graphene-encapsulated cobalt nanoparticles, which could possess efficient and stable electrocatalytic properties. Wei et al. [28] reported a reed-stalk-based Si-Fe/N/C ORR catalyst inspired by the Si-promoting graphitization findings. Deng et al. [29] described an activated carbon sheet derived from chitin, which was an efficient non-metal bifunctional electrocatalyst for both ORR and oxygen evolution reaction (OER). Li et al. [30] designed three-dimensional, honeycomb-like, nitrogen-doped carbon nanosheet/graphene nanonetwork films, which were prepared by dissolution and coagulation of chitin and graphene oxide in NaOH/urea aqueous solution using a repeated freezing–thawing process. Deng et al. [31] reported a chitin-derived porous carbon loaded with Co, N, and S heteroatoms, which showed excellent performance towards ORR, OER, and hydrogen evolution reaction (HER).

As numerous works reported, synthesizing non-precious ORR catalysts with high catalytic activity often requires further treatments such as activation [32] or hydrothermal [33] and chemical vaporous deposition [34]. However, intricate production steps are unpractical in industry and cannot achieve mass production. Herein, a one-step pyrolysis method was developed for preparation of a Co and N co-doped biocarbon catalyst derived from chitin. The obtained biomass, carbon-based material exhibits high catalytic activity and superior durability, which is comparable with the state-of-the-art noble-metal catalyst (Pt/C). The outstanding ORR performance can be attributed to the synergetic effect of the N-content functional group, the Co nanoparticles, and the carbon framework. Furthermore, an Al-air battery was assembled with a Co and N co-doped biocarbon electrocatalyst as an air cathode, which achieved high power density (32.24 mW cm−2) and excellent energy storage performance, indicating the Co and N co-doped biocarbon catalyst could easily be utilized as an efficient electrocatalyst in the Al-air battery.

#### **2. Results and Discussion**

#### *2.1. Morphological and Structural Characterization*

The morphology of as-prepared catalysts was first characterized by SEM and TEM. As shown in Figure 1a and Figure S1, a large number of pores are randomly distributed in the carbon skeleton of NC and CoNC. Furthermore, numerous particles are distributed in and on the skeleton of CoNC. In the TEM image shown in Figure 1b, several nanoparticles ranging from 15 nm to 40 nm can be observed in CoNC. From HR-TEM in Figure 1c, these nanoparticles are encapsulated in the CoNC skeleton and possess lattice spacing of 0.212 nm corresponding to Co particles. Raman spectra were applied to assess the defects of as-prepared materials. As shown in Figure 1d, Raman spectra of the NC and CoNC both have two obvious peaks (denoted as D band and G band) located in 1350 cm−1 and 1590 cm−1, which correspond to local defects of carbon materials and graphite degree graphitic in-plane vibration of the carbon materials, respectively. The intensity ratio of D band and G band (ID/IG) indicates the defect degree of carbon materials. As is shown in Figure 1g, the ID/IG of the CoNC and NC is about 1.06 and 0.95, respectively, which suggests that the CoNC materials contain more defects in the carbon skeleton than those of NC. Subsequently, the difference in morphology was further discussed via a nitrogen adsorption-desorption measurement. Based on nitrogen adsorption-desorption isotherms (Figure 1e), it can be noted that the CoNC possessed higher specific surface areas (~165 m2 g−1) than that of NC (~106 m2 g−1). Likewise, the CoNC also showed the lowest average pore size (2.99 nm), lower than that of NC (5.47 nm) (Figure 1f). The possible reason for the above phenomenon might be ascribed to

the sharp shrinking of the carbon matrix and catalytic effect induced by transition metal during the pyrolysis process [35,36].

**Figure 1.** The SEM and TEM images of cobalt nitrogen co-doped carbon material (CoNC) (**a**–**c**). Raman spectra (**d**), N2 adsorption/desorption isotherms (**e**) and the corresponding pore size (**f**).

The chemical states of as-prepared catalysts were studied by X-ray photoelectron spectroscopy (XPS). As shown in Figure 2a, the spectrum of CoNC contains four significant peaks with locations at 285 eV, 533 eV, 400 eV, and 782 eV corresponding to the C, O, N, Co elements, respectively. These results corresponded to the observation of the SEM and TEM images, suggesting that there are Co particles in CoNC. In a high-resolution C1s spectra of CoNC (Figure 2b), the detected peaks can be deconvoluted into 284.6 eV, 286.4 eV, 289.1 eV, and 285.4 eV, and each of them indicates the existence of C-C bond, C=O bond, O-C=O, and C-N bond in CoNC, respectively. Among these element peaks, the existence of a C-N bond implies the successful doping of N elements into the carbon skeleton. Subsequently, the types of N doped in carbon materials can be studied by high-resolution N1s spectra. As shown in Figure 2c, the N1s spectra can be fitted well with 398.4 eV, 399.2 eV, and 402.2 eV, which represent pyridinic-N, pyrrolic-N, and graphitic-N, respectively. It should be noted that the pyridinic-N is the first place where carbon-based ORR catalyst absorbs O2 and starts the ORR [37–40]. Meanwhile, the graphitic-N can become the electron-donor and facilitate charge mobility. As shown in Figure 2c and Figure S2, the CoNC material has a higher proportion of pyridinic-N and graphitic-N in N species than that of the NC material, which could further suggest that the CoNC is able to provide more efficient ORR performance than the NC.

**Figure 2.** X-ray photoelectron spectroscopy (XPS) survey spectrum of CoNC (**a**), C1s and N1s XPS spectrum of CoNC (**b**,**c**).

#### *2.2. Electrocatalytic Characteristics and Active Sites*

The electrocatalytic activities of the four as-prepared samples were first investigated via cyclic voltammetry (CV) tests performed via an electrochemical workstation (CHI 750e) with three-electrode systems. As shown in Figure 3a and Figure S3, NC and CoNC show obvious redox peaks located at 0.69 V and 0.81 V vs. RHE at O2-saturated KOH solution, and no obvious redox peak can be observed in CV curves at Ar-saturated KOH solution. This indicates that these redox peaks are formed by the ORR process. The CoNC catalyst possesses the positive ORR redox peak, which is comparable to the redox peak (0.82 V vs. RHE) shown in Pt/C (Figure 3b).

**Figure 3.** Cyclic voltammetry (CV) curves of CoNC (**a**) and Pt/C (**b**) in O2- and Ar-saturated 0.1 M KOH solution.

To further study the catalytic activities and ORR kinetics of the as-prepared catalysts, linear sweep voltammetry (LSV) tests were characterized with an rotating disk electrode (RDE) electrode. As shown in corresponding LSV curves (Figure 4a,b, the bar plot corresponds to current density and the line corresponds to onset potential), CoNC shows the higher positive onset potential (0.86V vs. RHE) and limiting current density (4.91 mA cm−2) than the PC (0.78 V, 4.33 mA cm−2) and NC (0.84 V, 4.32 mA cm−2) catalysts. Furthermore, these properties of CoNC are comparable to the benchmark catalyst of Pt/C, for which the onset potential and limiting current density are shown as 0.9 V and 5.67 mA cm−2, respectively.

As numerous works reported, ORR has two major pathways to reduce oxygen into water—the two-electron pathway and the four-electron pathway, respectively. In the four-electron pathway, oxygen is directly reduced to water and is able to provide more positive potential during the ORR process, whereas the two-electron pathway first reduces O2 to H2O2 then further reduces it to water. Thus, in order to investigate the catalytic pathway for CoNC, the ORR kinetics of CoNC were studied with RDE at different rotating speeds ranging from 100 rpm to 2500 rpm (Figure 4c). As shown in Figure 4d, LSV curves combined with the Koutecky-Lecich (K-L) equation are able to obtain several K-L plots at different potentials, which exhibits good parallelism and linearity, suggesting a first-order reaction toward ORR. The calculated electron transfer number is 3.73, which closed to the theoretical value 4, and suggests the four-electron dominant pathway during ORR. Meanwhile, Pt/C also shows the four-electron dominant pathway as electron transfer number (n = 4.03) calculated from K-L plots (Figure S4). For further study in ORR kinetic for CoNC, the rotating ring-disk electrode (RRDE) was used for testing the H2O2 productivity and electron transfer numbers while catalyzing ORR. As shown in Figure S5, CoNC produces a low proportion of H2O2 (less than 11.8%) and high electron transfer numbers (about 3.84) when ORR is carried out, while Pt/C shows 5.4% of H2O2 productivity and a 4.02 electron transfer number. These results elaborate that CoNC possesses the remarkable ORR catalytic activities in alkaline media and has the potential to be applied in metal-air batteries. The long-time durability tests were carried out by chronoamperometry measurement at 900 rpm in O2-saturated 0.l M KOH. As shown in Figure S6, after a 15000 s long-time durability test, the CoNC is

still able to provide 94.82% of initial current density, while the Pt/C can only maintain 76.64% of initial current density. This indicates that CoNC holds excellent long-term durability, which further benefits the stability of batteries assembled with CoNC.

**Figure 4.** The linear sweep voltammetry (LSV) curves of PC, NC, CoNC and Pt/C catalysts at 1600 rpm (**a**).The onset potentials and limiting current density of catalysts (**b**).The LSV curves of CoNC catalyst at the rotation speeds of 100, 225, 400, 625, 900, 1225, 1600, 2025 and 2500 rpm (**c**). The corresponding Koutecky-Levich (K-L) plots (**d**).

Subsequently, the home-built Al-air batteries were assembled with CoNC and Pt/C in order to demonstrate the practical application performance of catalysts. As shown in Figure 5a, the polarization curves of Al-air batteries were first used for estimating the capability of Al-air batteries. The corresponding power density calculated from the polarization curve of CoNC shows that the maximum power density it can provide is 32.24 mW cm−2, while the maximum power density for Al-air batteries assembled with Pt/C is 43.70 mW cm−2. The electrochemical impedance spectroscopy (EIS) was used to analyze the resistance in Al-air batteries. As calculated by the Nyquist plots in Figure S7, the internal resistance of the Al-air battery with CoNC is 3.46 Ω, which is comparable to that of the battery with 20 wt% Pt/C (3.29 Ω). In order to assess the long-term working property of the battery, the galvanostatic discharge performance was tested in a constant current of 20 mA cm−2 (Figure 5b). After a 5000 s test, the voltages of batteries with CoNC and Pt/C decreased 0.04 V and 0.21 V, respectively, and the voltages of batteries with CoNC remained at the constant discharge voltage of 1.17 V at 20 mA cm−2. The superior durability of batteries assembled with CoNC can be attributed to the excellent catalytic durability of CoNC.

**Figure 5.** The polarization curves and corresponding power densities of the Al-air batteries with two catalysts (**a**). The galvanostatic discharge curves of Al-air batteries with two catalysts at the constant current density of 20 mA cm−2 (**b**).

#### **3. Experimental Section**

#### *3.1. Preparation of Co and N Co-Doped Biocarbon*

The Co and N co-doped biocarbon is synthesized by a one-step pyrolysis method. Specifically, a mixture of chitin (2 g, MYM biological technology company limited, Nanjing, China), CoCl2·H2O (0.6 g, Fuchen chemical reagents factory, Tianjin, China) and urea (2 g, Beijing Chemical Works, Beijing, China) were ground thoroughly and dissolved into 40 mL of deionized water. Then the mixture was stirred overnight with a string speed of 600 rpm to form a homogeneous solution that was then dried at 60 ◦C for 48 h. Subsequently, the resulting sample was pyrolyzed at 800 ◦C for 1 h at a heating rate of 5 ◦C min−1 (Sager SG-GS1700). The resultant powder, N and Co co-doped biocarbon material, was denoted as CoNC. As controls, directly carbonized chitin (PC) and N-doped bio-carbon (NC) without any Co, were also prepared by the same method. For comparison, 20 wt% Pt/C was prepared and the loading for Pt was about 20 μg cm−2.

#### *3.2. Physicochemical Characterization*

The morphologies of as-prepared catalysts were studied by scanning electron microscopy (SEM) with a JEOL JSM-7500F field emission scanning electron microscope and transmission electron microscope (TEM) with a JEOL-2010 electron microscope operating at 200 kV. X-ray photoelectron spectroscopy (XPS) was performed on ESCALAB 250Xi (Thermo Scientific, Shanghai, China). A laser microscope system (B&W TEK BTC 162E) with an exciting laser of 532 nm was applied to study the Raman shifts. Then, the specific surface area and pore size were studied via automated gas sorption analyzer (ASAP 2020M-Physisoprtion Analyzer, Changchun, China) at 77 K. The Brunauer–Emmett–Teller method was used to calculate the specific surface area and the pore size.

#### *3.3. Electrochemical Measurements*

The electrochemical measurements were carried on a CHI 750E electrochemical workstation system with a typical three-electrode cell in 0.1 M KOH solution at room temperature. Glassy carbon electrode, Ag/AgCl, and Pt wire were used as a working electrode, reference electrode, and counter electrode, respectively. The catalyst inks were prepared by ultrasonically dispersing 5 mg of as-prepared catalysts, 5 mg of carbon black, and 50 μL of Nafion in 1 mL of ethanol. Then, 20 μL well-dispersed catalyst inks were loaded on a working electrode (5.5 mm diameter), leading to a 100 μg cm−2 catalyst loading. The as-prepared electrodes were dried completely at room temperature before the electrochemical tests. Cyclic voltammetry (CV) was measured in O2- and Ar-saturated 0.1 M KOH electrolyte with a scan rate of 10 mV s−1. For the ORR kinetics study, the linear sweep voltammetry (LSV) tests were

performed at rotating rates from 100 to 2500 in O2-saturated electrolyte with a scan speed of 10 mV s−1. The long-time durability of catalysts was tested in O2-saturated 0.l M KOH by chronoamperometry.

The electrons transfer number during ORR was analyzed from RDE measurements by using the following Koutecky-Lecich (K-L) equations (Equations (1) and (2)):

$$j^{-1} = j\_K^{-1} + j\_L^{-1} = j\_K^{-1} + B^{-1} \alpha^{-1/2} \tag{1}$$

$$\mathbf{B} = 0.2 \text{FC}\_0(\mathbf{D}\_0)^{2/3} v^{-1/6} \tag{2}$$

where *jL* and *jK* are the limiting current for the electrode reaction of reactive species by the diffusioncontrolled process and the kinetic current, respectively. *F* (96486.4 C mol−1) is the Faraday constant, n is the electron transfer number, C0\* (1.2 × 10−6 mol cm−3) is the concentration of O2 in 0.1 M KOH solution, D0 (1.9 × 10−5 cm2 s−1) is the diffusion coefficient of O2 in 0.1 M KOH, ω (rad s−1) is rotation rate and v (0.01 cm2 s−1) is the kinematic viscosity of the electrolyte [41,42].

#### *3.4. Al-Air Battery Tests*

The measurements of Al-air batteries were tested via home-built electrochemical cells. Al plate (99% purity) and 6 M KOH aqueous solution were used as the anode and electrolyte, respectively. A cathode was prepared by catalyst layer, carbon cloth (middle layer), carbon-based layer, and the diffusion layer. The mixture of carbon black and 40% Polytetrafluoroethylene (PTFE) was coated onto carbon cloth to form a carbon-based layer that prevents penetration of PTFE in the following process. A diffusion layer was prepared by applying a 60% PTFE layer coating on the carbon-based layer and heating it at 370 ◦C for 20 min. Subsequently, a catalyst layer was prepared by loading a mixture of CoNC, 5% Nafion, and ethanol on the other side of the carbon cloth (25 mg cm−2). For comparison, the Al-air battery that used commercial Pt/C as an air cathode was fabricated via the same procedures. The polarization curves and long-term discharge performances were tested by electrochemical workstation (CHI 750E) under air atmosphere.

#### **4. Conclusions**

In conclusion, a low-cost method was developed to synthesize efficient ORR catalyst derived from chitin. The biomass-derived electrocatalyst (CoNC) possesses high specific surface area, which contributes to enhanced ORR performance, and shows good ORR catalytic activities during ORR in alkaline media. In addition, Al-air batteries equipped with CoNC electrocatalyst show good power density and excellent long-term discharging durability, which is attributed to the efficient catalytic activities of CoNC. This work paves a new way for mass production of low-cost ORR catalyst and the utilization of biomass product.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4344/9/11/954/s1, Figure S1: The SEM and TEM images of NC, Figure S2: N1s XPS spectrum of NC, Figure S3: CV curves of NC catalysts, Figure S4: The LSV curves of Pt/C catalyst at the rotation speeds of 100, 225, 400, 625, 900, 1225, 1600, 2025 and 2500 rpm (a). The corresponding Koutecky-Levich (K-L) plots (b), Figure S5: The H2O2 yield and n value of CoNC and Pt/C, Figure S6: The i-t chronoamperometric curves of CoNC and Pt/C for 15000 s in an O2-saturated 0.1 M KOH solution, Figure S7: Nyquist plots of CoNC and commercial Pt/C from 100 KHz to 0.01 Hz.

**Author Contributions:** G.L. and S.Y. designed experiments; M.W. and J.M. carried out experiments; M.W. and Z.C. analyzed sequencing data; M.W. and H.Y. wrote the manuscript.

**Funding:** This work was funded by the National Natural Science Foundation of China (grant number: 51605188, 51605187), Department of Education of Jilin Province (grant number: JJKH20180093KJ, JJKH20181163KJ), and Jilin Provincial Science & Technology Department (Grant Number:<PHONE_NUMBER>4GX).

**Conflicts of Interest:** The author declares that there is no conflict of interests regarding the publication of this article.

#### **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*


**Chapter 9**

*New Advances on Fermentation Processes*

fungal microbiomes and mycotoxin production during ensiling and upon aerobic exposure. Applied and Environmental Microbiology.

[89] Dunière L, Xu S, Long J, et al. Bacterial and fungal core microbiomes associated with small grain silages during ensiling and aerobic spoilage. BMC Microbiology. 2017;**17**:50

[90] Ni K, Zhao J, Zhu B, et al. Assessing the fermentation quality and microbial community of the mixed silage of forage soybean with crop corn or sorghum. Bioresource Technology.

[91] Ren F, He R, Zhou X, et al. Dynamic changes in fermentation profiles and bacterial commjnity composition during sugarcane top silage fermentation: A preliminary study. Bioresource Technology. 2019;**285**:121315

2018;**84**:e02274-e02217

2018;**265**:563-567

**176**