**Optimisation, Optimal Control and Nonlinear Dynamics in Electrical Power, Energy Storage and Renewable Energy Systems**

Editors

**Victor M. Becerra Ahmed Rachid**

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

*Editors* Victor M. Becerra University of Portsmouth UK

Ahmed Rachid University of Picardie Jules Verne France

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

This is a reprint of articles from the Topic published online in the open access journals *Energies* (ISSN 1996-1073), *Electronics* (ISSN 2079-9292), *Batteries* (ISSN 2313-0105), and *Electricity* (ISSN 2673-4826) (available at: https://www.mdpi.com/topics/optimisation).

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**, *Volume Number*, Page Range.

**ISBN 978-3-0365-5585-0 (Hbk) ISBN 978-3-0365-5586-7 (PDF)**

Cover image courtesy of Victor M. Becerra

© 2022 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 Editors**

### **Victor M. Becerra**

Professor Victor M. Becerra is currently a Professor of Power Systems Engineering at the School of Energy and Electronic Engineering, University of Portsmouth, UK, where he has been since 2015. He obtained his PhD in Control Engineering from the City University of London in 1994. During 1989–1991, he carried out power systems analysis and control systems design at C.V.G. Edelca, Caracas, Venezuela. Between 1994 and 1999, he was a post-doctoral Research Fellow at the Control Engineering Research Centre at the City University of London. During the period between 2000 and 2015, he was an academic at the School of Systems Engineering, University of Reading, UK, where he became a Professor of Automatic Control in 2012. His current research interests consider a range of issues related to power systems and smart grids, as well as the methods and applications of automatic control. During his career, he has received funding for research and innovation from the EPSRC, Innovate UK, the Royal Academy of Engineering, the European Space Agency, and the EU. He has published over 150 refereed papers, one research monograph and two edited volumes. Professor Becerra is a Fellow of the IET and a Senior Member of the IEEE.

#### **Ahmed Rachid**

Ahmed Rachid received an engineering degree in Cybernetics in 1983, an MBA in 1985, a PhD in electrical engineering in 1986 and a Habilitation degree in 1991 from the University of Nancy (France). He has been teaching since 1984 in different universities and engineering schools mainly in the fields of control engineering and signal processing. Since 1992, he has been a professor at the University of Picardie Jules Verne (France), where he created the institution's first laboratory in Automatic Systems. His current research includes processes and systems modeling, simulation, control, diagnosis and observation with applications in blast furnaces, vehicle dynamics, diesel engines, induction electrical drives, stepper motors, greenhouses, magnetic levitation, pantograph-catenary, co-generation, electric vehicles, solar power systems, etc. He has supervised 24 PhD theses and has co-authored 2 patents, 2 books and over 150 scientific papers. He has organized and chaired several international conferences and has coordinated various European research projects and industrial R&D contracts.

## *Article* **Hardware Implementation of Novel Shade Dispersion PV Reconfiguration Technique to Enhance Maximum Power under Partial Shading Conditions**

**Vinaya Chandrakant Chavan 1 , Suresh Mikkili 1,\* and Tomonobu Senjyu 2**

**Abstract:** Partial shade condition is a significant factor contributing to the PV panel performance in mismatch losses and power generation. The technique suggested in this study allows the physical rearrangement of the PV panel to distribute the shade on the entire PV array. MPPT, selecting suitable inverter topology, or PV panel reconfiguration enhances the performance of the PV panel. This study proposes a new shade dispersing method, novel shade dispersion (NSD). It compares the performance of the NSD method with conventional configurations (CCs). This research article models and simulates 6 × 6 PV array configurations such as Series-Parallel (SP), Total-Cross-Tide (TCT), Bridge-Linked (BL), Honey-Comb (HC), and the newly proposed NSD method under nonshading and nine different partial shading cases. The performance indices used for comparative analysis are global maximum power points, efficiency, power enhancement, open circuit voltage, short circuit current, and number of crests. The Soltech 1 STH-215-P PV module was selected in the MATLAB/Simulink environment to simulate PV array arrangements. Hardware experiments validate the performance of the NSD method.

**Keywords:** partial shading; conventional configuration; reconfiguration techniques

## **1. Introduction**

Renewable energy has become one of the most critical factors to preserve the climate and the earth's assets for the coming generation. Power has a crucial part in the financial growth and the welfare of a nation. From India's perspective, solar energy increases power generation and generates energy reliably considering ecological, communal, and financially beneficial properties. India has significant potential, availability, and other features of solar power. Therefore, the government also emphasizes the promotion of photovoltaic generation in the Indian power sector [1]. Series and parallel connected modules form an array to obtain power per requirement. PV material, temperature, irradiance, dust, array structure, maximum power point tracking (MPPT) method, converter topology, etc., are the different parameters responsible for the performance of PV array [2,3]. This uncertainty affects the estimation of a PV cell's series and shunt resistances (namely, Rs and Rsh). Ref. [4] provided the uncertainty analysis of the measured currents, voltages, temperatures, and irradiances. At a constant temperature, the output current increases with an increase in irradiance with negligible effect on the voltage. Voltage is substantially reduced at constant irradiance, and the current slightly increases with increased temperature up to the rated temperature [5]. Any object that obstructs the sun rays' path to the PV panel or uneven falling sunlight on the earth creates a partial shading condition (PSC) [6]. A bypass diode is connected antiparallel to the PV module to avoid the hot spot issues created by PSC. Under PSC, the difference between the minimum irradiance received by non-overlapped

**Citation:** Chavan, V.C.; Mikkili, S.; Senjyu, T. Hardware Implementation of Novel Shade Dispersion PV Reconfiguration Technique to Enhance Maximum Power under Partial Shading Conditions. *Energies* **2022**, *15*, 3515. https://doi.org/ 10.3390/en15103515

Academic Editor: Alberto Reatti

Received: 7 April 2022 Accepted: 9 May 2022 Published: 11 May 2022

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

**Copyright:** © 2022 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/).

and overlapped solar cells determines the overcurrent [7]. A diode adds antiparallel either to a PV module or a set of PV series modules [8]. It inserts the number of peaks in the PV curve, among which the rise consisting of the highest value is known as the global maximum power point (GMPP). Traditional MPPT struggles to obtain true MPP in rapidly changing PSCs. In this area, extensive research considers improved conventional methods, new metaheuristic methods, etc., which increases the complexity and cost of the MPPT system [9]. Ref. [10] performed the experimental validation of TCT with the standard deviation of PSC and found a better configuration. Results demonstrated that 15% and 25% reduced irradiance. The power yield degraded to 4% and 7.5%. Refs. [11,12] explored the performance of series (S), parallel (P), SP, TCT, BL, and HC PV topologies and concluded, without considering bypass diode BL, and with bypass diode TCT generates maximum power. Reconfiguration is an effective way to improve power generation [13]. This was classified into two broad categories: electrical and physical reconfiguration. Both methods have their boons and downsides [14]. Ref. [15] proposed the maximum and minimum algorithm dynamic reconfiguration method, which uses a double pole double throw switch to change the electrical connections. Ref. [16] provided the dynamic array reconfiguration for water pumping applications. The authors of [17] investigated efficiency and robustness of the dynamic array reconfiguration for hybrid PV system application. Refs. [15–17] efficiently used the switch matrix to reconfigure the PV modules in the array for different applications. Ref. [18] proposed an electrical repositioning on a 4 × 5 arrangement with 12 switches to make adaptive cross ties to enhance the power generation. For case I the proposed method improved GMPP, fill factor (FF), and power enhancement (PE) compared with SP, and minimized power losses (PL). It was found best at 291.3 W, 0.586, 2.71%, 55.6 W, and 83.97%, for TCT in MATLAB/Simulink and experimentation with the same performance index under all PSCs. Complexity, sensor requirement, switches, and their controlling structure, sensor speed dependency, cost, etc., are the significant limitations of dynamic reconfiguration. Therefore, many researchers go for static array reconfiguration. The authors of [19] proposed and demonstrated a Sudoku-based physical reconfiguration system on a 9 × 9 array. It improves the power generation from 4 to 26% for different PSCs. The authors of [20] proposed an improved Sudoku method on a 9 × 9 array and compared it to the CC. Ref. [21] offered a magic square method and compared it with TCT, Sudoku, and the optimal Sudoku method on an 81 modules array structure. Ref. [22] proposed the zigzag method and compared it with the TCT and Sudoku method on a 6 × 6 PV array. Ref. [23] related competence square and the dominance square method with TCT on a 9 × 9 PV array. Ref. [24] proposed a clue-based skyscraper puzzle method and compared it with TCT, Sudoku method, and the puzzle-based competence square method. This study reviews many static reconfiguration (SR) methods and analyzes their performance. SR is a robust, simple, less costly, and effective way to deal with partial shading issues and improve the performance of the PV array. Still, these methods have some limitations.

The novelty of work is:

Electrical networks persist similarly to TCT and the physical position of units in the PV array changed.

• The proposed scheme was analyzed under ten shading cases. The results validated with SP, TCT, BL, and HC compared the global maximum power point (GMPP), mislead power (MP), output power (%O/P), efficiency, fill factor, and % power enhancement for NSD under-considered shading pattern. Figure 1 gives the workflow of this paper. same column, maintaining the least possible wiring; • The proposed scheme was analyzed under ten shading cases. The results validated with SP, TCT, BL, and HC compared the global maximum power point (GMPP), mislead power (MP), output power (%O/P), efficiency, fill factor, and % power enhancement for NSD under-considered shading pattern. Figure 1 gives the workflow of this paper.

ray. This method considers magic sub-matrix shifting to reposition the panels in the

*Energies* **2022**, *15*, x FOR PEER REVIEW 3 of 18

**Figure 1.** Representation of PV array connection structures under partial shading conditions with performance evaluation parameters. **Figure 1.** Representation of PV array connection structures under partial shading conditions with performance evaluation parameters.

#### **2. Modeling and Simulation of NSD Method 2. Modeling and Simulation of NSD Method**

This section explains the modeling and simulation of conventional PV array configurations and the proposed novel shade dispersion (NSD) method in MATLAB/Simulink. The panel electrical parameters has mentioned in Appendix A. The flowchart for panel repositioning is provided in Figure 2. The present work consists of 36 modules placed in six rows and six column patterns (shown in Figure 3a) for conventional configurations and the NSD method (Figure 3b). The TCT connection-based MATLAB simulation is shown in Figure 4. Mathematical modeling for output voltage, output current, and output power of all the traditional PV array structures are provided in Table 1. Here, the module, row, and array output voltages are V୨, Vୖ, and V, respectively. The module, string, and array output currents are I୨, Iୗ, and I. P is the PV generation of the array. A novel shade dispersion method is proposed to reconfigure the panel of a variety to improve the performance through shade dispersion. In the proposed technique, electrical connections remain the same as in TCT [20]; only the physical location of the modules changes. It uses voltage, current, and power equations for TCT, as This section explains the modeling and simulation of conventional PV array configurations and the proposed novel shade dispersion (NSD) method in MATLAB/Simulink. The panel electrical parameters has mentioned in Appendix A. The flowchart for panel repositioning is provided in Figure 2. The present work consists of 36 modules placed in six rows and six column patterns (shown in Figure 3a) for conventional configurations and the NSD method (Figure 3b). The TCT connection-based MATLAB simulation is shown in Figure 4. Mathematical modeling for output voltage, output current, and output power of all the traditional PV array structures are provided in Table 1. Here, the module, row, and array output voltages are Vj , VR, and VPV, respectively. The module, string, and array output currents are Ij , IS, and IPV. PPV is the PV generation of the array. A novel shade dispersion method is proposed to reconfigure the panel of a variety to improve the performance through shade dispersion. In the proposed technique, electrical connections remain the same as in TCT [20]; only the physical location of the modules changes. It uses voltage, current, and power equations for TCT, as mentioned in Table 1. The basic principle used for this method is row current equalization. It proposes a magic submatrix

*Energies* **2022**, *15*, x FOR PEER REVIEW 4 of 18

Step I: Formation of the submatrix

Step II: Formation of the primary matrix

the same row and vacant column serially;

the magic number 21;

fifth step is achieved.

is 6;

mentioned in Table 1. The basic principle used for this method is row current equalization. It proposes a magic submatrix arrangement with the magic number 21 and its logi-

1. The submatrix is proposed so that rows and column multiplication of the submatrix

2. Therefore, numbers 1 to 6 are placed in the 2 × 3 submatrix to summation equal to

3. The submatrix placing is started in the 1st row, 1st column position. The next submatrix is placed at the next (3rd) row, next (2nd) column position. This process con-

4. Submatrix rows are interchanged, and that submatrix starts to fill in vacant places of the primary matrix, starting from the next column of previously filled areas. Here, the interchanged submatrix is filled from 1st row, 4th column. The next submatrix is placed in the 3rd row and 4th column. The last column of the submatrix is placed in

cal shifting to reposition the panels within a column, as shown in Figure 3b. The steps involved to build the proposed 6 × 6 NSD matrix are:

tinues until the last (6th) row is filled by the considered submatrix;

arrangement with the magic number 21 and its logical shifting to reposition the panels within a column, as shown in Figure 3b. 5. It is performed so that the summation of row numbers, as well as column numbers, is equal to the magic number 21, as shown in Figure 3b; 6. The process continues until vacant places in the primary matrix are filled, and the

**Figure 2.** Flowchart for panel repositioning of NSD method (**a**) Step 1: submatrix formation; (**b**) Step 2: main matrix formation. **Figure 2.** Flowchart for panel repositioning of NSD method (**a**) Step 1: submatrix formation; (**b**) Step 2: main matrix formation.

C1 C2 C3 C4 C5 C6 R1 **111111** R2 **222222** R3 **333333** R4 **444444** R5 **555555** R6 **666666**

> **264 531**

Submatrix

**21**

(a)

(b)

**Figure 3.** Panel positions in PV array (**a**) conventional methods, (**b**) NSD method.

**11 12 13 14 15 16 21 22 23 24 25 26 31 32 33 34 35 36 41 42 43 44 45 46 51 52 53 54 55 56 61 62 63 64 65 66**

C1 C2 C3 C4 C5 C6 R1 **264531** 21 R2 **531264** 21 R3 **126453** 21 R4 **453126** 21 R5 **312645** 21 R6 **645312** 21 21 21 21 21 21 21

Step 2: main matrix formation.

(**b**) **Figure 2.** Flowchart for panel repositioning of NSD method (**a**) Step 1: submatrix formation; (**b**)

**Stop**

**Yes**

**Are all the places filled?**

**Yes Yes**

**No No**

**Start**

**Select square matrix of a×a**

**Start placing submatrix from 1st row, 1st column**

> **Is 'a' prime number**

> > **Place submatrix of m×n from next row, next column of previous submatrix**

**No**

**Is 'ath ' row filled**

**Place next submatrix element rowwise staring from next row and jump by skipping one column of previous submatrix**

**Yes**

**No**

**Is 'ath ' row filled**

**Figure 3.** Panel positions in PV array (**a**) conventional methods, (**b**) NSD method. **Figure 3.** Panel positions in PV array (**a**) conventional methods, (**b**) NSD method.

**Figure 4.** TCT-based MATLAB simulation for the NSD method. **Figure 4.** TCT-based MATLAB simulation for the NSD method.

**Table 1.** The output voltage, current, and power for conventional PV array configurations. **Table 1.** The output voltage, current, and power for conventional PV array configurations.

**HC** V = V୨ ୨ୀଵ = 6V୨ I = Iୗଵ Iୗଶ Iୗଷ Iୗସ Iୗହ Iୗ = 6Iୗ P = 36V୨Iୗ The steps involved to build the proposed 6 × 6 NSD matrix are: Step I: Formation of the submatrix

global maximum power point (GMPP), mislead nine PSCs as power (MP), output power (%O/P), efficiency, fill factor, and % power enhancement of NSD. The %PE was calcu-

• Case II to case IV have group 1 of 24 panels with 1000 W irradiance; groups 2, 3, and 4 have 200 W, 500 W, and 700 W irradiance, respectively, with four panels in each

• In case V, 10 modules have PSCs, group 1 with 1000 W has 26 panels. Group 2, with 200 W, has four modules, and groups 3 and 4 have 500 W, 700 W irradiance, respec-

• Case VI has groups 1, 2, 3, and 4 with 1000 W, 200 W, 500 W, and 700 W irradiance

• In case VII, 24 panels came in group 1. Groups 2 and 3, have 200 W and 700 W irradiance, respectively. Group 4 has two panels that have 700 W irradiance; • In case VIII, groups 1, 2, 3, and 4 have five panels in each group of 1000 W, 200 W,

500 W, and 700 W irradiance with 16, 8, 7, and five panels, respectively;

non-shading condition), the performance of the NSD method is the same as the CC (i.e., without reconfiguration), as shown in Figure 7a. The performance indicators are the Step II: Formation of the primary matrix

on 28, 3, 1, and 4 modules, respectively;

group;

tively, with three panels;

#### **3. Partial Shadings Conditions and Shade Dispersion Effect**

For CCs the shade does not scatter, as shown in Figure 5. Through NSD, the partial shading condition spread on a complete array, as shown in Figure 6. Under case I (i.e., non-shading condition), the performance of the NSD method is the same as the CC (i.e., without reconfiguration), as shown in Figure 7a. The performance indicators are the global maximum power point (GMPP), mislead nine PSCs as power (MP), output power (%O/P), efficiency, fill factor, and % power enhancement of NSD. The %PE was calculated for NSD and compared with all the CCs. For case I, all 36 modules of the PV array have 1000 W irradiance; the remaining PSCs have four groups of irradiance. *Energies* **2022**, *15*, x FOR PEER REVIEW 7 of 18 • In case IX, only five modules are partially shaded, and 31 modules have 1000 W irradiance in group 1. It has two 200 W modules in group 2, one 500 W module in group 3 and two 700 W irradiance panels in group 4; • In case X, with 1000 W, group 1 has 25 panels. Group 2, with 200 W has four panels, and groups 3 and 4 have 500 W and 700 W irradiance, respectively, with three panels. *Energies* **2022**, *15*, x FOR PEER REVIEW 7 of 18 • In case IX, only five modules are partially shaded, and 31 modules have 1000 W irradiance in group 1. It has two 200 W modules in group 2, one 500 W module in group 3 and two 700 W irradiance panels in group 4; • In case X, with 1000 W, group 1 has 25 panels. Group 2, with 200 W has four panels, and groups 3 and 4 have 500 W and 700 W irradiance, respectively, with three pan-

**Figure 5.** Partial shading conditions (**a**) case II, (**b**) case III, (**c**)case IV, (**d**) case V, (**e**) case VI, (**f**) case VII, (**g**) case VIII, (**h**) case IX, (**i**) case X. **Figure 5.** Partial shading conditions (**a**) case II, (**b**) case III, (**c**)case IV, (**d**) case V, (**e**) case VI, (**f**) case VII, (**g**) case VIII, (**h**) case IX, (**i**) case X. (g) (h) (i) **Figure 5.** Partial shading conditions (**a**) case II, (**b**) case III, (**c**)case IV, (**d**) case V, (**e**) case VI, (**f**) case VII, (**g**) case VIII, (**h**) case IX, (**i**) case X.

**1000 1000 1000<PHONE_NUMBER>**

**1000<PHONE_NUMBER> 1000 1000**

**1000 700 1000 1000 200 1000**

VI, (**f**) case VII, (**g**)case VIII, (**h**) case IX, (**i**) case X. **Figure 6.** Shade spreading due to NSD method (**a**) case II, (**b**) case III, (**c**)case IV, (**d**) case V, (**e)** case VI, (**f**) case VII, (**g**)case VIII, (**h**) case IX, (**i**) case X. **Figure 6.** Shade spreading due to NSD method (**a**) case II, (**b**) case III, (**c**)case IV, (**d**) case V, (**e)** case VI, (**f**) case VII, (**g**)case VIII, (**h**) case IX, (**i**) case X.

**Figure 7.** PV curves under various irradiance levels (**a**) case I, (**b**) case II, (**c**)case III, (**d**) case IV, (**e**) case V, (**f**) case VI, (**g**)case VII, (**h**) case VIII, (**i**) case IX, (**j**) case X. **Figure 7.** PV curves under various irradiance levels (**a**) case I, (**b**) case II, (**c**)case III, (**d**) case IV, (**e**) case V, (**f**) case VI, (**g**)case VII, (**h**) case VIII, (**i**) case IX, (**j**) case X.

#### **4. Result and Discussion**

The performance evaluation of PV array CCs and proposed reconfiguration was performed under various shading situations on MATLAB/Simulink. The assessment has two segments, (1) power-voltage profile and (2) computational parameters.

#### *4.1. Power-Voltage Profiles under Various PSC*

The sketches in Figure 7 represent the variation in PV curves for different conventional configurations and the proposed NSD method of reconfiguration under various irradiance levels. Under case I of the non-shading condition, all the panels have the same standard irradiance of 1000 W. Figure 7a highlights that all the configuration and NSD methods have the same ideal performance under case I with only one peak on the PV curve.

Under case II, conventional configurations have two peaks on the PV curve, while the NSD method has only one extreme. TCT, BL, and HC have identical extreme power, and it is more than SP, while for NSD, the maximum capacity is much more enhanced, as revealed in Figure 7b. Figure 7c indicates that SP has the lowest power output for case III. It offers continuous improvement for BL, HC, TCT, and NSD. NSD offers 6234 W as the maximum power output. In Figure 7d, SP has three peaks, and other methods have two heights. BL and HC show improved performance over SP but lesser than TCT and NSD. NSD provides the highest output of 6234 W, as shown. Figure 7e portrays that BL and HC have six peaks, SP and TCT have five peaks, while NSD is reduced to four peaks and the highest power generation equal to 6303 W for case V. In Figure 7f case VI, HC and SP have the same maximum power generation as BL and TCT with three peaks; in contrast, the NSD has the highest power generation of 6646 W with three small deep peaks. Figure 7g showcases the case VII PV curves for SP, TCT, BL, HC, and NSD. NSD has a single top with 6133 W as the highest power generation. In Figure 7h SP and TCT have six and four peaks. BL, HC, and NSD structures have three peaks with the highest power, 5160 W, for the NSD method. Figure 7i portrays case IX. Under this shade, TCT has maximum power extraction of 6563 W in conventional configuration while NSD has improved power extraction with 6953 W.

#### *4.2. Computational Parameters*

The performances of all CC and NSD were investigated in the following six computational parameters, as observed in Appendix B, which highlights the efficacy of NSD over CC. Tables 2–6 give the performance of SP, TCT, BL, HC and NSD method under ten PSCs respectively.

**Table 2.** SP performance under PSC.

**Table 3.** TCT performance under PSC.

**Table 4.** BL performance under PSC.

**Table 5.** HC performance under PSC.

**Table 6.** NSD performance under PSC.

### 4.2.1. Global Maximum Power Point (GMPP)

The PV curve has multiple peaks due to PSC, as shown in Figure 8. GMPP is the maximum power of all the crests. Figure 8a indicates that GMPP for all CCs and NSD have concentrated 7665 W. For case II, TCT, BL, and HC provide equal power at 5039 W in CC, and NSD improves it up to 6234 W. Case III also NSD has GMPP at 6234 W. Under case IV in CC, TCT provides maximum output at 6028 W while NSD has 6211 W. Under case V, NSD has 6303 W as a GMPP; this is 1201 W more than HC, which offers the highest power of all CCs. Under case VI to case X in CC, TCT offers the uppermost GMPP while NSD has GMPP at 6646 W, 6133 W, 5160 W, 6953 W, and 6350 W for respective PSC. It is 812 W, 1325 W, 988 W, 389 W, and 257 W more than TCT under the same PSCs.

#### 4.2.2. Mislead Power (MP)

Figure 8b shows that the proposed NSD method provides the least misleading power as the red NSD trace is at the innermost on the radar pattern. Case I does not have any MP. From case II to case X it is 1431 W, 1431.8 W, 1454 W, 1363 W, 1020 W, 1533 W, 2506 W, 713 W, and 1315 W continuously. NSD reduces the MP by about 11 to 46% under these shading cases.

#### 4.2.3. Percent Output (%O/P)

Figure 8c Highlights that all considered structure provides 100% output for case I. NSD provides 81% O/P, the maximum from cases II, III, and IV. In CC, HC provides the highest %O/P at 66.5% under case V, while NSD improves it to 82.69%. In CC under case VI to case X, TCT offers the uppermost %O/P, and NSD improved it by 10.5%, 17.3%, 12.9%, 5%, and 3.4% compared with TCT for the respective shading cases.

#### 4.2.4. Efficiency (η) −

Figure 8d shows the NSD structure is more efficient than TCT, which has the same electrical connections; under-considered PSCs. The radar pattern representation in Figure 8d clearly shows that NSD offers efficiency near 21%. It is approximately equal to the PV array efficiency for case I.

#### 4.2.5. Fill Factor (FF)−

Figure 8e highlights that the proposed method offers 2–5% improvement in FF for cases III, IV, IX, and X; 15.19% more FF for case V, 8% improvement for case VI, and 20–25% enhancement for case II, VII, and VIII, compared with maximum FFs in CC.

#### 4.2.6. Percent Power Enhancement (%PE)

The %PE ensures the power increment by the proposed reconfiguration. Figure 8f shows that NSD improves power extraction compared with CC. Compared with SP, TCT, BL, and HC, NSD has a 25.9, 15.79, 17.44, and 16.45% power boost.

4.2.5. Fill Factor (FF)−

4.2.6. Percent Power Enhancement (%PE)

Figure 8e highlights that the proposed method offers 2–5% improvement in FF for cases III, IV, IX, and X; 15.19% more FF for case V, 8% improvement for case VI, and 20–

The %PE ensures the power increment by the proposed reconfiguration. Figure 8f

25% enhancement for case II, VII, and VIII, compared with maximum FFs in CC.

BL, and HC, NSD has a 25.9, 15.79, 17.44, and 16.45% power boost.

**Figure 8.** Performance evaluation parameter (**a**) GMPP (W), (**b**) MP(W), (**c**)%O/P, (**d**) η, (**e**) FF, (**f**) **Figure 8.** Performance evaluation parameter (**a**) GMPP (W), (**b**) MP(W), (**c**)%O/P, (**d**) η, (**e**) FF, (**f**) %PE.

#### **5. Experimental Implementation and Result Analysis 5. Experimental Implementation and Result Analysis**

#### *5.1. Experimental Setup 5.1. Experimental Setup*

%PE.

Towards real-time implementation of the proposed NSD method, a prototype laboratory experiment with 16 PV modules and a 4 × 4 array structure was conducted on the PV chroma 62050 H-600 S simulator (Taiwan). The chroma 63203 E DC (Taiwan) electronic load was used to develop an experimental setup, as shown in Figure 9. The electrical parameters for each 5 Watt panel are, VOC = 21.6 V, VM = 17.5 V, ISC = 0.3 A, IM = 0.29 V. Under non-shading conditions, this PV array provides 80 W power generation. Towards real-time implementation of the proposed NSD method, a prototype laboratory experiment with 16 PV modules and a 4 × 4 array structure was conducted on the PV chroma 62050 H-600 S simulator (Taiwan). The chroma 63203 E DC (Taiwan) electronic load was used to develop an experimental setup, as shown in Figure 9. The electrical parameters for each 5 Watt panel are, VOC = 21.6 V, VM = 17.5 V, ISC = 0.3 A, IM = 0.29 V. Under non-shading conditions, this PV array provides 80 W power generation.

> As shown in Figures 10 and 11, two types of irradiance circumstances were examined in the experiment. The shading patterns considered for experimentation were horizontal (case I) and triangular corner (case II). These are similar to case II and case V in Section 2 but with different irradiance levels and array sizes. The shade concentrates before reconfiguring this array structure, as shown in Figure 10, for all the conventional configurations. After applying the proposed NSD reconfiguration, the shadow was distributed on the array, as shown in Figure 11.

**Figure 9.** Experimental setup. **Figure 9.** Experimental setup. configurations. After applying the proposed NSD reconfiguration, the shadow was distributed on the array, as shown in Figure 11. configurations. After applying the proposed NSD reconfiguration, the shadow was distributed on the array, as shown in Figure 11.

*Energies* **2022**, *15*, x FOR PEER REVIEW 13 of 18

*Energies* **2022**, *15*, x FOR PEER REVIEW 13 of 18

**1000 1000 1000 1000 1000 1000 1000 1000 Figure 10.** PSC before reconfiguration (**a**) case I, (**b**) case II. **Figure 10.** PSC before reconfiguration (**a**) case I, (**b**) case II. **Figure 10.** PSC before reconfiguration (**a**) case I, (**b**) case II.

**200 1000 1000 1000 500 1000 700 1000 Figure 11.** PSC after reconfiguration (**a**) case I, (**b**) case II. **Figure 11.** PSC after reconfiguration (**a**) case I, (**b**) case II. **Figure 11.** PSC after reconfiguration (**a**) case I, (**b**) case II.

#### **1000 300 1000 1000** *5.2. Experimental Result Analysis 5.2. Experimental Result Analysis 5.2. Experimental Result Analysis*

**1000 1000 1000 500 1000 1000 400 1000<PHONE_NUMBER> 1000 300 1000 1000 1000** a) b) To perform the analysis of conventional configurations, SP, BL, HC, TCT, and the proposed NSD reconfiguration, the PV curve nature was obtained by changing the load resistance through DC Electronic load from 4000 Ω to 1 Ω. The curve nature was observed on an HDO4034A (USA) 350 MHz high-definition oscilloscope. To perform the analysis of conventional configurations, SP, BL, HC, TCT, and the proposed NSD reconfiguration, the PV curve nature was obtained by changing the load resistance through DC Electronic load from 4000 Ω to 1 Ω. The curve nature was observed on an HDO4034A (USA) 350 MHz high-definition oscilloscope. To perform the analysis of conventional configurations, SP, BL, HC, TCT, and the proposed NSD reconfiguration, the PV curve nature was obtained by changing the load resistance through DC Electronic load from 4000 Ω to 1 Ω. The curve nature was observed on an HDO4034A (USA) 350 MHz high-definition oscilloscope.

**1000 700 1000 1000**

#### **Figure 11.** PSC after reconfiguration (**a**) case I, (**b**) case II. 5.2.1. PSC Case I

*5.2. Experimental Result Analysis*  To perform the analysis of conventional configurations, SP, BL, HC, TCT, and the proposed NSD reconfiguration, the PV curve nature was obtained by changing the load resistance through DC Electronic load from 4000 Ω to 1 Ω. The curve nature was observed Figure 12 depicts the PV curves for each conventional configuration, NSD reconfiguration scheme, and the global MPP (the dot on the PV curve) under PSC case I. The PV curve shows GMPP for SP, BL, HC, TCT, and NSD at 59.22 W, 49.79 W, 49.79 W, 59.23 W, and 66.57 W, respectively. Moreover, SP and TCT have two peaks on the PV curve, BL and HC have three heights, and NSD reconfiguration reduces the number of peaks.

on an HDO4034A (USA) 350 MHz high-definition oscilloscope.

Figure 12 depicts the PV curves for each conventional configuration, NSD reconfiguration scheme, and the global MPP (the dot on the PV curve) under PSC case I. The PV curve shows GMPP for SP, BL, HC, TCT, and NSD at 59.22 W, 49.79 W, 49.79 W, 59.23 W, and 66.57 W, respectively. Moreover, SP and TCT have two peaks on the PV curve, BL and HC have three heights, and NSD reconfiguration reduces the number of peaks.

**Figure 12.** PV curves under PSC case I. **Figure 12.** PV curves under PSC case I.

#### 5.2.2. PSC Case II 5.2.2. PSC Case II

5.2.1. PSC Case I

Under this shading condition, Figure 13 shows that all conventional configurations have three peak PV curves. In contrast, the proposed NSD reconfiguration has one minor local peak and GMPP at a higher level with an improved curve nature. SP, BL, HC, TCT, and NSD provide GMPP at 52.92 W, 54.63 W, 53.75 W, 56.04 W, and 70.03 W, respectively. Figure 13 highlights that there is much improvement in the PV curve nature under this shading condition and GMPP location. Under this shading condition, Figure 13 shows that all conventional configurations have three peak PV curves. In contrast, the proposed NSD reconfiguration has one minor local peak and GMPP at a higher level with an improved curve nature. SP, BL, HC, TCT, and NSD provide GMPP at 52.92 W, 54.63 W, 53.75 W, 56.04 W, and 70.03 W, respectively. Figure 13 highlights that there is much improvement in the PV curve nature under this shading condition and GMPP location. *Energies* **2022**, *15*, x FOR PEER REVIEW 15 of 18

**Figure 13.** PV curves under PSC case II. **Figure 13.** PV curves under PSC case II.

#### **6. Conclusions 6. Conclusions**

spectively;

From the discussion mentioned above, this study concludes that: From the discussion mentioned above, this study concludes that:

• In partial shading cases, for the least irradiating case VIII and most irradiating case IX, NSD has maximum power extraction of 5160 W and 6952.8 W, respectively; the

• Under PSC for the considered PV array, TCT has the highest GMPP, %O/P, efficiency, FF, %PE, and least %MP at 5683 W, 74%, 18.54%, 0.5784, and 28.73%. Compared with TCT, NSD has 705.77 W, 9.2%, and 0.1139 enhancement in GMPP, %O/P,

• The results obtained from the first simulation show that the NSD method achieved the most significant output power and overall performance compared with all con-

• Moreover, experimental research was carried out. The suggested NSD PV array reconfiguration design beats the other conventional PV array configurations (SP, BL, HC, and TCT) in power increase and several peaks. The experimental setup vali-

**Author Contributions:** Conceptualization—V.C.C. and S.M.; Methodology, Software, Validation—V.C.C., S.M. and T.S.; Investigation, Resources, Writing—Original Draft Preparation, Writing—Review and Editing, Visualization; Supervision; Project Administration, V.C.C., S.M. and T.S.

**Funding:** The authors gratefully acknowledge the support offered by the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India, under the

dated the influence of panel shifting on various PV array topologies.

ventional configurations under nine partial shading cases;

All authors have read and agreed to the published version of the manuscript.

and FF, and a 2.41% reduction in MP;

Grant number: EEQ/2021/000294 for this research work. **Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable. **Data Availability Statement:** Not applicable.

**Author Contributions:** Conceptualization—V.C.C. and S.M.; Methodology, Software, Validation— V.C.C., S.M. and T.S.; Investigation, Resources, Writing—Original Draft Preparation, Writing—Review and Editing, Visualization; Supervision; Project Administration, V.C.C., S.M. and T.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors gratefully acknowledge the support offered by the Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India, under the Grant number: EEQ/2021/000294 for this research work.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **Appendix A**

**Table A1.** PV Module electrical characteristics.

#### **Appendix B**

**Definitions**:

• **Global Maximum Power Point (GMPP):** GMPP is the single maximum point on the PV curve.

#### GMPP = max(multiple power Peaks)

• **Mislead Power (MP)**: MP denotes the difference between GMPP in the non-shading (NS) condition and GMPP in the shaded (PSC). It creates the hot spot effect and, due to that, power loss. It is provided by:

MP = GMPPNS − GMPPPSC

• **Percent Output (%O/P):** The %O/P is calculated as the maximum power ratio under PSC to extreme power under non-shading (NS) conditions. It provides how much energy is generated under PSC compared with power at standard irradiance. It reaches the %O/P of the PV array under PSC and NS conditions. It provides how much percentage of PSC affects the PV array output.

$$\% \text{O/P} = \frac{\text{P}\_{\text{max}(\text{PSC})}}{\text{P}\_{\text{max}(\text{NS})}}$$

• **Efficiency (**η**):** It is provided by the ratio output power (power generation) POut to input power (solar insolation) Pin On the PV panel. It is provided by:

$$\eta = \frac{\mathcal{P}\_{\text{Out}}}{\mathcal{P}\_{\text{in}}}$$

• **Fill Factor (FF):** The maximum power generated at PSC to the total rated capacity of a solar PV array is FF. It is also known as the usefulness of PV array. It is always 0 to 1. It is provided by:

$$\text{FF} = \frac{\text{V}\_{\text{M}}\text{I}\_{\text{M}}}{\text{V}\_{\text{OC}}\text{Is}}$$

• **Percent Power Enhancement (%PE)**: %PE is provided by the difference between the GMPP for the proposed NSD method and the GMPP for conventional configuration (CC) under partial shading conditions. It provides how much power is enhanced by NSD compared with CC under various irradiance conditions.

$$\% \text{PE} = \frac{\text{GMPP}\_{\text{NSD}} - \text{GMPP}\_{\text{CC}}}{\text{GMPP}\_{\text{CC}}}$$

#### **References**

## *Article* **Multiobjective Load Dispatch for Coal-Fired Power Plants under Renewable-Energy Accommodation Based on a Nondominated-Sorting Grey Wolf Optimizer Algorithm**

**Yue Cao 1,\*, Tao Li 1 , Tianyu He 1 , Yuwei Wei 2 , Ming Li 3 and Fengqi Si 1**

**Abstract:** Coal-fired power plants are widely used to achieve a power balance in grids with renewable energy, which leads to new requirements for speediness in load dispatch. This paper presents a nondominated-sorting grey wolf optimizer algorithm (NSGWO) for the multiobjective load dispatch of coal-fired power plants that employed efficient nondominated sorting, a reference-point selection strategy, and a simulated binary crossover operator. The optimization results of the benchmark functions indicated that the NSGWO algorithm had a better accuracy and a better distribution than the traditional multiobjective grey wolf optimizer algorithm. Regarding the load dispatch of economy, environmental protection, and speediness strategies, the NSGWO had the best performance of all the simulated algorithms. The optimal-compromise solutions of the economy and speediness strategies of the NSGWO algorithm had a good distribution, which elucidated that this novel algorithm was favorable to allowing coal-fired power plants to accommodate renewable energy.

**Keywords:** load dispatch; multiobjective optimization; grey wolf optimizer; nondominated sorting; coal-fired power plant

## **1. Introduction**

Modern coal-fired power plants are clean and highly efficient, and their emissions of NOx and coal consumption are controlled at remarkably low levels [1]. Regarding the national energy distribution in China, coal-fired power plants still need to output power for basic loads in the grid. A coal-fired power plant may have different units (such as water-cooling and air-cooling units [2]) and units of various scales. Differences in unit characteristics may result in variable performances; therefore, it is necessary to optimize the operation strategy of coal-fired power plants with multiple units on the basis of both the economy and environmental-protection objectives. Furthermore, regarding the longterm objective of carbon neutrality [3], coal-fired power plants will play a new role in the power grid. Considering the development of renewable energy including solar power and wind power, coal-fired power plants, which have a strong capability to accommodate renewable energy, are widely used to achieve a power balance in power grids. In terms of the uncertainty of renewable energy [4], it is necessary to operate coal-fired power plants under conditions of rapid load change. Therefore, the optimal operation of coal-fired power plants should be further investigated while considering the factors of the economy, environmental protection, and speediness.

**Citation:** Cao, Y.; Li, T.; He, T.; Wei, Y.; Li, M.; Si, F. Multiobjective Load Dispatch for Coal-Fired Power Plants under Renewable-Energy Accommodation Based on a Nondominated-Sorting Grey Wolf Optimizer Algorithm. *Energies* **2022**, *15*, 2915. https://doi.org/10.3390/ en15082915

Academic Editor: Elyas Rakhshani

Received: 21 March 2022 Accepted: 13 April 2022 Published: 15 April 2022

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

**Copyright:** © 2022 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/).

#### *1.1. Load-Dispatch Strategy*

Units of various types and scales in a coal-fired power plant usually receive the load command from the grid, which limits the optimal operation of the coal-fired power plant. Regarding the load-dispatch strategy, the coal-fired power plant receives the total load command from the grid, and then dispatches it to different units [5]. Based on previously mentioned factors, namely the economy, environmental protection, and speediness, the load dispatch of a coal-fired power plant is a multiobjective optimization of load command between different units.

Several researchers have employed optimal load-dispatch strategies for power generation systems. Dai et al. [6] proposed a multiobjective economic load-dispatch method for large coal-fired power plants using a hierarchical clustering and retrieval strategy based on the fuzzy c-means clustering algorithm. They concluded that the data-mining-based method was able to achieve plant-level optimal load dispatching while meeting the actual requirements of the grid. Ishraque et al. [7] introduced various optimization strategies of load dispatch for microgrid power systems with renewable energy sources. They considered that the proposed load-following method was the best dispatch strategy to achieve lower operating costs and lower pollutant emissions. Cui et al. [8] considered that the heat load and power load could be optimized simultaneously in microgrids. Regarding the differences in real-time response of heat load and power load, they performed short-term and long-term simulations and elucidated that the total cost could be reduced, as could the peak load of the power grid. Li et al. [9] also focused on investigating multiobjective load dispatch between the heat load and power load, and concluded that the proposed method could balance the production cost and pollutant-emission objectives under the fluctuation of the thermoelectric load demand. Xu et al. [10] introduced the load-dispatch strategy to save coal for coal-fired power plants using a support-vector machine. Their results indicated that the proposed load-dispatch strategy performed better than the normal separate load demand on a single unit in a coal-fired power plant. In addition, novel loaddispatch strategies have been studied for similar power-generation systems. Liao et al. [11] investigated the short-term load-dispatching method for a hydropower plant with multiple turbines, and simulation results indicated that it could achieve optimal operation while simultaneously meeting the practical operation requirements of the hydraulic and electric constraints. Jamal et al. [12] introduced a bioinspired computational heuristic algorithm for dispatch loads for a power-generation system that included coal-fired thermal-generating units and wind-power units. They considered that the proposed algorithm had advantages in accuracy, convergence, and robustness in solving load-dispatch problems. Furthermore, Bie et al. [13] conducted a load-dispatch investigation based on the uncertainty of renewable energy, and the results showed that the real-time load-dispatch model could calculate the load command while considering the primary and secondary regulation of the grid.

In general, much research has been published on load-dispatch optimization, but few studies considered the present application of a coal-fired power plant with renewableenergy accommodation, which is significant in terms of the multiobjective load dispatch of modern coal-fired power plants.

#### *1.2. Multiobjective Optimization Algorithms*

Among the load-dispatch strategies, researchers consider the optimization algorithm to be the key technology needed to achieve the objectives of energy saving and reducing pollutant emissions. In addition to the traditional optimization methods (including dynamic programming, Lagrange relaxation, and the exhaustive method), intelligent optimization algorithms have already been introduced in recent investigations.

Intelligent optimization algorithms usually employ an evolutionary mechanism or heuristic algorithm to enhance the convergence speed and robustness ability of the optimal solution procedure. Ghasemi et al. [14] proposed a differential evolution algorithm (DE) based on the mathematical model of sociopolitical evolution. They performed comparative simulations to verify the effectiveness of the algorithm, and concluded that this algorithm

was a reliable tool for solving the load-dispatch problem. Neto et al. [15] proposed a method that combined a differential-evolution algorithm and the greedy randomized adaptive search procedure algorithm. The simulation results showed that the algorithm could solve the nonsmooth load-dispatch problem by considering the effect of valve-point loading. Secui et al. [16] proposed a new heuristic algorithm called modified symbiosis (MSOS) that could solve the large-scale-unit load-optimization problem with a valve-point effect. They concluded that the proposed algorithm had a good performance based on the simulation results of five power-generation systems of different constraints and dimensions.

Furthermore, other heuristic algorithms, including the bat algorithm and the population extremum algorithm, have been introduced to solve the multiobjective optimization of load dispatch. Kavousi-Fard et al. [17] proposed an improved bat algorithm to solve the nonconvex load-optimization problem, which was then examined on IEEE test systems to show its high abilities and satisfying performance. Chen et al. [18] used the feasible-solution method to deal with constraints, and proposed a constrained multiobjective population extremal optimization algorithm called CMOPEO-EED. They considered the proposed algorithm to be a viable option to solve the load-dispatch problem with renewable power generation.

Particle-swarm optimization (PSO) is another commonly used algorithm in multiobjective optimization. Hosseinnezhad et al. [19] used an improved PSO algorithm (θ-PSO) to solve the economic load-optimization problem that could easily deal with various constraints. The findings of this study elucidated that the advantage of computational efficiency of the θ-PSO algorithm could be a promising alternative approach to solving the load-dispatch problem in practical power systems. Regarding the combined advantages of the PSO algorithm and other intelligent algorithms, researchers have investigated hybrid algorithms based on PSO. Narimani et al. [20] proposed a hybrid algorithm based on the PSO algorithm and leapfrog algorithm for multiregional economic and environmentalprotection load distribution. They verified the effectiveness of the hybrid algorithm using different test systems to balance the multiobjectives of generation costs and environmental issues. Zou et al. [21] proposed a new global particle-swarm-optimization algorithm, which adopted an updated method that relied on the global optimal particle to guide all particle-search activities and a slightly perturbed randomization method for particle flight trajectories with uniform distribution. Mandal et al. [22] introduced variable acceleration coefficients into the PSO algorithm and then proposed a self-organized hierarchical particleswarm-optimization algorithm. They found that the results of the proposed method were superior in terms of fuel costs, emission output, and losses.

In general, intelligent optimization algorithms have advantages in solving the multiobjective optimization problems of load dispatch at power plants, though their application under the conditions of renewable energy accommodation should be further investigated.

#### *1.3. Grey Wolf Optimizer*

Former studies show that various intelligent algorithms have been employed to solve the multioptimization issues in power-generation systems and similar applications. However, considering the rapid load-change conditions present under renewable-energy accommodation, the load dispatch of current coal-fired power plants requires rapid and accurate convergence.

The grey wolf optimizer (GWO) is a swarm-intelligence optimization algorithm proposed by Mirjalili et al. [23]. Compared with other natural heuristic optimization algorithms such as particle-swarm optimization and differential evolution, GWO has attracted attention due to its fast convergence speed, strong robustness, and high solution accuracy with a simple structure. Zhang et al. [24] used a GWO to plan the path of an unmanned combat aerial vehicle, and they concluded that it was more competent than other state-of-the-art evolutionary algorithms in terms of quality, speed, and stability of solutions. Although the grey wolf algorithm has superior performance on most data sets, it may not have the best performance on all. In view of the defects of a GWO, many scholars have made

relevant improvements to the performance of the GWO algorithm. Madhiarasan et al. [25] optimized the initial population of a GWO to improve the diversity of the population. Saremi et al. [26] improved the search mechanism of a GWO. The results indicated that the proposed GWO algorithm was able to significantly improve the performance of multiobjective optimization. Long et al. [27] proposed an updating method for the key parameter of a GWO in the application of constrained optimization problems. They considered that the simulations for both benchmark functions and engineering applications could demonstrate the good performance of the proposed GWO algorithm. Moreover, hybrid optimization algorithms combining GWO with other intelligent optimization algorithms have already been introduced by researchers [28–30]. However, few scholars have conducted in-depth research on the evolution mechanism of the GWO for the present application of loaddispatch optimization. Song et al. [31] introduced a GWO to solve the economic emission problem of two different power systems, and used a price penalty factor to combine the goals of economic dispatch and emission dispatch into a single function. Wong et al. [32] solved the economic scheduling problem under practical constraints with a GWO and tested them on two test systems with practical constraints. Moradi et al. [33] used a GWO to solve nonlinear and nonconvex economic scheduling problems while considering the valve-point effect and transfer loss, and compared the obtained results with some existing heuristic methods. Jangir et al. [34] investigated the economic constrained-emission dispatch problem based on a nondominated-sorting GWO, and they considered that it could be used in a power grid with wind power. However, they did not consider the speediness requirement of units in the coal-fired power plant, which indicated that the solving procedure of nondominated-sorting GWO should be further optimized.

The above shows that limited work has been conducted on the multiobjective optimization algorithm while considering the rapid and accurate convergence requirements under the unit rapid load-change condition. The more constraints that are considered, the more complex the load-dispatch optimization will be. Therefore, further research on the GWO needs to be carried out according to the specific problems of the multiobjective optimization of the load dispatch of coal-fired power plants.

#### *1.4. Present Investigation*

On the basis of the previous literature review, although much research has been published on load-dispatch optimization, few studies considered the present application of a coal-fired power plant with a renewable-energy accommodation. Moreover, limited work has been conducted on the multiobjective optimization algorithm that considered the rapid and accurate convergence requirements under the unit's rapid load-change condition. Therefore, an investigation into the multiobjective optimization of the load dispatch of a coal-fired power plant was conducted using a novel grey wolf algorithm.

In the present paper, the proposed nondominated-sorting grey wolf algorithm showed that it may help solve various optimization problems in similar applications. The organization of this study was arranged as follows. Section 2 shows the methods for the multiobjective optimization problem of load dispatch. Section 3 details the solving procedure of the novel grey wolf algorithm. Section 4 presents the simulation results and further discussion. Section 5 presents the conclusions of this study.

#### **2. Methods**

A theoretical model for the multiobjective optimization of the load dispatch of a coalfired power plant, as well as optimization constraints and performance evaluation indices, were proposed in order to achieve the requirements related to the economy, environmental protection, and speediness.

#### *2.1. Multiobjective Optimization Problem of a Coal-Fired Power Plant*

#### 2.1.1. Problem Description

The load-dispatch problem of a coal-fired power plant can generally be described as per Equation (1):

$$\begin{cases} \min \mathbf{F}(\mathbf{x}) = \min \left( \mathbf{f}\_1(\mathbf{x}), \mathbf{f}\_2(\mathbf{x}), \dots, \mathbf{f}\_m(\mathbf{x}) \right)^T\\ \text{s.t.} \begin{cases} \mathbf{g}\_i(\mathbf{x}) \le \mathbf{0}, i = 1, 2, \dots, p\\ \mathbf{h}\_j(\mathbf{x}) = \mathbf{0}, j = 1, 2, \dots, q \end{cases} \end{cases} \tag{1}$$

where F(*x*) is the objective function of the economy performance, and the sub-functions represent the economy characteristics of units in the coal-fired power plant. Further, g*i* (*x*) and h*j*(*x*) represent inequality constraints and equality constraints, respectively. Specifically, pollutant emission is an inequality constraint when it is not selected as an optimization objective.

For a coal-fired power plant, the relationship between the pollutant emission g*i* and the unit load *Pi* can be expressed as the quadratic function relationship shown below

$$\mathbf{g}\_i(P\_i) = \mathbf{a}\_i P\_i^2 + \beta\_i P\_i + \gamma\_i \tag{2}$$

where *αi* , *βi ,* and *γi* are the emission-characteristic coefficients of the unit.

Regarding the requirements of the relevant national policies of environmental protection, it is necessary to employ the environmental protection index as an optimization object instead of an inequality constraint. Thus, the objective function of load dispatch when considering environmental protection is defined as shown in Equation (3):

$$\min \mathcal{G} = \min \sum\_{i=1}^{N} \mathbf{g}\_i(P\_i) \tag{3}$$

where G represents the pollutant discharge of the whole plant.

The power grid has specific requirements for the load response rate of coal-fired power plants, specially under the condition of renewable-energy accommodation. Therefore, the speediness of load variation should be considered in the load dispatch of coal-fired power plants. The minimum time consumption of load dispatch under ideal conditions *Tideal* is shown in Equation (4):

$$T\_{ideal} = \left| P\_D - \sum\_{i=1}^{N} P\_{i,now} \right| / \sum\_{i=1}^{N} v\_i \tag{4}$$

where *PD* is the total load command of all units, MW; and *vi* represents the load changing speed of unit *i*, MW/min.

The sum of squares of the differences between the load-dispatch times of each unit and the minimum time under ideal conditions T(*P*) can be calculated as follows:

$$\text{T}(P) = \sum\_{i=1}^{N} \left( (P\_i - P\_{i,now}) / v\_i - T\_{ideal} \right)^2 \tag{5}$$

Therefore, combining former indicators with the economic indicators, the multiobjective optimization of load dispatch for a coal-fired power plant can be described as shown in Equation (6):

$$\begin{array}{l} \min(F, G, T), \min(F, G), \min(F, T) \\ \text{s.t.} \begin{cases} \sum\_{i=1}^{N} P\_i - P\_D \\ P\_i^{\min} \le P\_i \le P\_i^{\max} \\ \left| P\_i - P\_i^{\max} \right| / t\_i \le v\_i \end{cases} \end{array} \tag{6}$$

#### 2.1.2. Constraint-Violation Degree

The constraint-violation degree is introduced to show whether an individual *x* violates the *i*-th inequality constraint or the *j*-th equality constraint, which is expressed in Equations (7) and (8):

$$\text{Con}\_i(\mathbf{x}) = \max\{\mathbf{g}\_i(\mathbf{x}), \mathbf{0}\}, 1 \le i \le p \tag{7}$$

$$\mathcal{Con}\_j(\mathbf{x}) = \max\{h\_j(\mathbf{x}), 0\}, 1 \le j \le q \tag{8}$$

Thus, the total constraint-violation degree of individual *x* is shown in Equation (9):

$$\nu(\mathbf{x}) = \sum\_{i=1}^{p} \mathbb{C}on\_i(\mathbf{x}) + \sum\_{j=1}^{q} \mathbb{C}on\_j(\mathbf{x}) \tag{9}$$

In the multiobjective optimization of the load dispatch for a coal-fired power plant, the equation constraint describes that the sum of the loads of each unit must be equal to the total load command, while the inequality constraint shows the load-regulation range.

#### 2.1.3. Feasible and Infeasible Solutions

When an individual *x* satisfies both the inequality constraints and the equality constraints in the load-dispatch model, this individual is called a feasible solution; otherwise, it is an infeasible solution. The region composed of all feasible solutions is called the feasible region Ω, and the set composed of infeasible solutions is called the infeasible region.

#### 2.1.4. Dominance

For the two specific solutions *x*1 and *x*2, the economy, speediness, and environmentalprotection indicators corresponding to solution *x*1 are better than those corresponding to solution *x*2. Thus, solution *x*1 is called a nondominated solution, while solution *x2* is called a dominant solution, which can be written as *x*1 *x*2.

All individuals in the solution set that cannot be dominated by other individuals are called nondominated solutions, and this set of nondominated solutions is called the Pareto front (PF) of optimization solutions.

#### *2.2. Constraints for Multiobjective Optimization*

Constraint processing should be conducted to balance the relationship between the objective function and the constraints. In this paper, the penalty-function method was used to deal with constraints. By adding the penalty function *F*(*x*, *M*) into the objective function, the constrained optimization problem was transformed into an unconstrained optimization problem. Parameter *M* in the penalty function represents a large positive number, called the penalty factor, which acts as a penalty for the solution that violates the constraints. Using the external penalty function, the multiobjective function with the penalty function is expressed as follows:

$$\min F(\mathbf{x}, M) = \min \left( f\_1(\mathbf{x}) + M\_1 v(\mathbf{x}), f\_2(\mathbf{x}) + M\_2 v(\mathbf{x}), \dots, f\_m(\mathbf{x}) + M\_m v(\mathbf{x}) \right) \tag{10}$$

#### *2.3. Evaluation Index of Multiobjective Optimization*

For the multiobjective optimization of a load dispatch, evaluation indices were employed to compare the simulation results, which may have shown the convergence performance and comprehensive performance of the optimization algorithm.

#### 2.3.1. Convergence Evaluation Index

The generational distance (*GD*) was introduced to evaluate the convergence performance, which represented the degree of the solution set *P* to the real Pareto frontier *P\**, as shown in Equation (11):

$$GD(P, P^\*) = \frac{1}{|P|} \sqrt{\sum\_{i=1}^{|P|} d\_i^2} \tag{11}$$

where |*P*| is the number of solutions in *P*, and *di* represents the Euclidean distance from the *i*-th solution in *P* to the nearest solution in *P\**. Therefore, a small *GD* value meant the convergence performance was better.

#### 2.3.2. Comprehensive Evaluation Index

The inverse generational distance (*IGD*) was used to represent the degree of the real Pareto frontier *P\** to the solution set *P* [35]. It evaluated the convergence and diversity of the algorithm simultaneously, which was defined as:

$$IGD(P, P^\*) = \frac{1}{|P^\*|} \sqrt{\sum\_{i=1}^{|P^\*|} d\_i'^2} \tag{12}$$

where *d* 0 *i* 2 represents the nearest Euclidean distance from the *i*-th individual in *P\** to the solution in *P*. Here, a small *IGD* showed that the optimization algorithm was better in both the convergence performance and the diversity of solutions.

#### 2.3.3. Optimal-Compromise Solution

For the multiobjective optimization of a load dispatch, the fuzzy set theory is suggested to find an optimal-compromise solution from the Pareto fronts. The fuzzy membership degree for the load-dispatch optimization can be expressed as:

$$\varphi\_{pi} = \begin{cases} 1 & f\_{pi} \le f\_i^{\min} \\ (f\_i^{\max} - f\_{pi}) / (f\_i^{\max} - f\_i^{\min}) & f\_i^{\min} < f\_{pi} < f\_i^{\max} \\ 0 & f\_{pi} \ge f\_i^{\max} \end{cases} \tag{13}$$

where *ϕpi* represents the satisfaction of the *i*-th objective of the *p*-th optimal solution, *i* ∈ {1, 2, . . . , *m*}; *fpi* represents the value of the *i*-th objective of the *p*-th optimal solution; and *f* max *i* and *f* min *i* represent the maximum and minimum values of the *i*-th objective in all optimal solutions, respectively. When *ϕpi* was equal to 1, it meant that the individual was completely satisfied with the *i*-th target value.

The standardized satisfaction calculation method of each solution *p* in the nondominated solution set is shown in Equation (14), where the solution corresponding to the largest *ϕpi* was the optimal-compromise solution:

$$\varphi\_p = \sum\_{k=1}^{m} \varphi\_{pk} / \sum\_{p=1}^{N} \sum\_{k=1}^{m} \varphi\_{pk} \tag{14}$$

Regarding the previous objectives and constraints of the multiobjective optimization of the load dispatch for a coal-fired power plant, the sequence diagram is shown in Figure 1. Figure 1.

Regarding the previous objectives and constraints of the multiobjective optimization of the load dispatch for a coal-fired power plant, the sequence diagram is shown in

**Figure 1.** Sequence diagram of multiobjective optimization of load dispatch of coal-fired power **Figure 1.** Sequence diagram of multiobjective optimization of load dispatch of coal-fired power plant.

#### plant. **3. Solving Procedure of the Optimization Algorithm**

**3. Solving Procedure of the Optimization Algorithm**  The multiobjective grey wolf algorithm was employed to find the load-dispatchoptimization solution for a coal-fired power plant. However, the issues of local optimum and insufficient stability may limit the performance of this optimization algorithm. In addition, the influence of the alpha wolf may result in a poor distribution of the obtained solutions., and the traditional nondominated sorting may consume more time in the iteration process. Considering the multiobjective optimization of a load dispatch with multiple complex constrains, a nondominated-sorting grey wolf optimizer (NSGWO) algorithm was proposed that employed efficient nondominated sorting (ENS), a reference-The multiobjective grey wolf algorithm was employed to find the load-dispatchoptimization solution for a coal-fired power plant. However, the issues of local optimum and insufficient stability may limit the performance of this optimization algorithm. In addition, the influence of the alpha wolf may result in a poor distribution of the obtained solutions., and the traditional nondominated sorting may consume more time in the iteration process. Considering the multiobjective optimization of a load dispatch with multiple complex constrains, a nondominated-sorting grey wolf optimizer (NSGWO) algorithm was proposed that employed efficient nondominated sorting (ENS), a referencepoint selection strategy, and a simulated binary crossover operator.

#### point selection strategy, and a simulated binary crossover operator. *3.1. Fast Nondominated Sorting Based on ENS*

lowing equation:

*3.1. Fast Nondominated Sorting Based on ENS*  Fast nondominated sorting is a cyclic stratification process. If this process is conducted in sequence, it will require too much time to converge. In fact, the frontier num-Fast nondominated sorting is a cyclic stratification process. If this process is conducted in sequence, it will require too much time to converge. In fact, the frontier number of individuals is determined by the dominated individual, which could infer the following equation:

$$front(p) = 1 + \max\_{q \in Q} front(q) \tag{15}$$

q Q *front front* (p) 1 max (q) ∈ = + (15) where *f ront*(*p*) represents the number of the frontier where individual *p* is located, and *Q* is the set composed of all individuals that dominate individual *p*.

where *front p*( ) represents the number of the frontier where individual *p* is located, and *Q* is the set composed of all individuals that dominate individual *p*. The first step of the ENS was to arrange all individuals in the population according to their first-dimension target values from small to large. If the first-dimension target values were the same, we compared their second-dimension index values, and so on. In the sorted population, it was impossible for the individual in the latter order to domi-The first step of the ENS was to arrange all individuals in the population according to their first-dimension target values from small to large. If the first-dimension target values were the same, we compared their second-dimension index values, and so on. In the sorted population, it was impossible for the individual in the latter order to dominate the individual in front of it. This meant that the individual that dominated *p* must have been assigned to the frontier when examining *p*. The number of the frontier where *p* is located was calculated using Equation (15).

nate the individual in front of it. This meant that the individual that dominated *p* must have been assigned to the frontier when examining *p*. The number of the frontier where *p* is located was calculated using Equation (15). On this basis, starting from the first frontier *F*1, it was judged in turn whether each frontier contained an individual that dominated *p*. The first frontier that does not contain the dominant *p* was the frontier where *p* resides. If there were individuals dominating p up to the last frontier *Fk* , then *p* was assigned to a new frontier *Fk*+1. If all individuals were mutually nondominant, then a nondominant comparison was required between them. In this case, since there was only one frontier, the complexity was *Tworst* = *m* · *N*(*N* −

1)/2 = *O*(*mN*2 ). The algorithm could achieve the best performance if the population satisfied the following conditions: there were √ *N* frontiers in total, and each frontier contained √ *N* individuals. Moreover, an individual was dominated by all individuals whose frontier number was smaller than its under this condition. For a specific individual *p*, if it belonged to frontier *Fi* , it was necessary to find an individual dominating it in each frontier *Fi*(*j* < *i*). This indicated that at least *i* − 1 nondominated comparisons were required. For all √ *N* individuals in *Fi* , a total of √ *N*(*i* + 1) comparisons were required. In addition, all individuals of √ *Fi* needed to be compared pairwise, which meant that a total of *N*( √ *N* − 1)/2 comparisons were required. Considering all frontiers, the best period of complexity is shown in Equation (16):

$$T\_{\rm best} = m \cdot \sum\_{i=1}^{\sqrt{N}} \left[ \sqrt{N}(i-1) + \sqrt{N}(\sqrt{N}-1)/2 \right] = m \cdot \left[ N(\sqrt{N}-1)/2 + N(\sqrt{N}-1)/2 \right] = O(\text{MN}^{1.5}) \tag{16}$$

#### *3.2. Selection Strategy Based on Reference Point*

To select *N* individuals from the set of offspring and parents R*t* = *Pt*∪*Qt* , it should firstly divide R*t* into multiple nondominated layers (*F*1, *F*2, . . . , *Fk* ) by nondominated sorting. Then, a new population *St* is established from the start of *F1* until its size is equal to or exceeds *N* for the first time. In order to maintain the diversity of solutions, a selection strategy based on reference point was used to select *nL* = *N* − ∑ *L*−1 *i*=1 *Ni* individuals, where *Ni* is the number of individuals in the *i*-th frontier, and *nL* is the number of individuals that needed to be selected from the last frontier satisfying the former requirements. Moreover, the systematic approach was introduced to generate a set of uniform weight vectors {*λ*1, *λ*2, . . . , *λR*} so that the reference points could be determined. In order to compare the target values of each dimension, it was necessary to adaptively normalize the target value of each dimension. The ideal point − z = (*z* min 1 , *z* min 2 , . . . , *z* min *m* ) of defining population *St* was composed of the minimum value *z* min 1 , *i* = 1, 2, . . . , *m* of each dimension of *St* . Based on the ideal point, the objective function could be transformed to *f* 0 *i* (*x*) = *fi*(*x*) − *z* min *i* . Then, a hyperplane could be built by extra points obtained from the objective function; this process can be expressed as:

$$\text{ASF}(\mathbf{x}, w) = \max\_{i=1}^{m} f\_i'(\mathbf{x}) / w\_{i\nu} \mathbf{x} \in \mathbb{S}\_t \tag{17}$$

$$Z^{i, \text{max}} = s: \text{argmin}\_{\text{s} \in \text{S}\_{\text{I}}} ASE(s, w), w = (\tau, \dots, \tau), \tau = 10^{-6} \tag{18}$$

For the *i*-th target, an extra target vector *z i*,max was generated. Thus, *m* extra vectors were generated by *m* targets which further constituted an *m*-dimensional linear hyperplane, and then the intercept *ai* , *i* = 1, 2, . . . , *m* could be calculated. The objective function can be normalized as shown in Equation (19):

$$f\_i^n(\mathbf{x}) = (f\_i'(\mathbf{x}) - z\_i^{\min}) / (a\_i - z\_i^{\min}) \tag{19}$$

#### *3.3. Simulated Binary Crossover Operator*

Under the influence of the top three levels of wolves, a GWO easily falls into a local optimum. Therefore, a simulated binary crossover (SBX) operator was introduced in this paper. Assuming that the two parent individuals are *x* 1 (*x* 1 1 , *x* 1 2 , . . . *x* 1 *n* ) and *x* 2 (*x* 2 1 , *x* 2 2 , . . . *x* 2 *n* ), the SBX operator could generate the offspring individuals *c* 1 (*c* 1 1 , *c* 1 2 , . . . *c* 1 *n* ) and *c* 2 (*c* 2 1 , *c* 2 2 , . . . *c* 2 *n* ) using the following equations:

$$\begin{cases} \ c\_i^1 = 0.5 \times \left[ (1+\beta) \cdot x\_i^1 + (1-\beta) \cdot x\_i^2 \right] \\\ c\_i^2 = 0.5 \times \left[ (1-\beta) \cdot x\_i^1 + (1+\beta) \cdot x\_i^2 \right] \end{cases} \tag{20}$$

where *β* is randomly determined by the distribution factor *η* according to Equation (21):

$$\beta = \begin{cases} \left(rand \times 2\right)^{1/(1+\eta)}, rand \le 0.5\\ \left(1/\left(2 - rand \times 2\right)\right)^{1/(1+\eta)}, rand > 0.5 \end{cases} \tag{21}$$

With the increase in iterations, the probability of generating individuals according to the grey wolf hunting mechanism increased, thereby enhancing the development ability and convergence of the algorithm. In this paper, the variable a was controlled by a nonlinear decreasing strategy to balance the exploration and development capabilities of the algorithm, as shown in Equation (22): With the increase in iterations, the probability of generating individuals according to the grey wolf hunting mechanism increased, thereby enhancing the development ability and convergence of the algorithm. In this paper, the variable a was controlled by a nonlinear decreasing strategy to balance the exploration and development capabilities of the algorithm, as shown in Equation (22):

Under the influence of the top three levels of wolves, a GWO easily falls into a local optimum. Therefore, a simulated binary crossover (SBX) operator was introduced in this

1 2 ( , ,... ) *n x xx x* , the SBX operator could generate the offspring individuals

1 12 2 12 0.5 [(1 ) (1 ) ] 0.5 [(1 ) (1 ) ] *i ii i ii*

 β

1/(1 )

+

η

(1/ (2 2)) , 0.5

*rand rand*

 β

(20)

*c xx c xx* β

β

where *β* is randomly determined by the distribution factor *η* according to Equation (21): 1/(1 )

( 2) , 0.5

+

η

−× >

*rand rand*

× ≤ =

1 2 ( , ,... ) *n ccc c* using the following equations:

 = × + +− = × − ++

paper. Assuming that the two parent individuals are 111 1

$$a = 2(1 - (t/iter)^2) \tag{22}$$

1 2 ( , ,... ) *n x xx x* and

(21)

where *iter* is the maximum number of iterations. where *iter* is the maximum number of iterations.

β

*Energies* **2022**, *15*, x FOR PEER REVIEW 10 of 20

*3.3. Simulated Binary Crossover Operator*

1 2 ( , ,... ) *n ccc c* and 222 2

222 2

111 1

#### *3.4. Solving Procedure of NSGWO 3.4. Solving Procedure of NSGWO*

The solving procedure of the proposed nondominated-sorting grey wolf optimizer algorithm is presented in Figure 2. There were five specific steps of this algorithm, which are detailed as follows. The solving procedure of the proposed nondominated-sorting grey wolf optimizer algorithm is presented in Figure 2. There were five specific steps of this algorithm, which are detailed as follows.

**Figure 2.** Solving procedure of nondominated-sorting grey wolf optimizer algorithm. **Figure 2.** Solving procedure of nondominated-sorting grey wolf optimizer algorithm.

Step 1: Set the number *N* of individuals in the grey wolf population and the maximum number of iterations *iter*. Initialize parameters *a*, *A,* and *C*. Step 1: Set the number *N* of individuals in the grey wolf population and the maximum number of iterations *iter*. Initialize parameters *a*, *A*, and *C*.

Step 2: Initialize the grey wolf population *P*, and generate a reference point *Z*.

Step 2: Initialize the grey wolf population *P*, and generate a reference point *Z*. Step 3: Calculate the target value of each dimension of individuals in the popula-Step 3: Calculate the target value of each dimension of individuals in the population.

tion. Step 4: Perform the nondominated sorting based on ENS, and select the first three levels of wolves *α*, *β*, *δ*.

Step 5: Update the population *Pt* and generate the offspring *Qt* . When the random number is greater than the set probability threshold, binary crossover mutation is conducted. Otherwise, update the population based on the evolution mechanism of GWO.

Step 6: Select the first *N* individuals based on the reference point in the merged population R*t* = *Pt*∪*Qt* , and then generate a new population *St* .

Step 7: If the maximum number of iterations is not reached, repeat steps (2–4). Otherwise, output the optimization results.

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

#### *4.1. NSGWO Algorithm Performance Test*

In order to test and evaluate the performance of the NSGWO algorithm, three benchmark functions [36] were employed in this paper, the detailed structures of which are shown in Table 1.

**Table 1.** Benchmark functions for multiobjective optimization test of NSGWO algorithm.

Then, test simulations were conducted for the multiobjective optimization algorithm; the results for the *GD* and *IGD* values of the former benchmark functions are presented in Tables 2 and 3, respectively. Here, the results of the multiobjective grey wolf optimizer (MOGWO) and the multiobjective particle-swarm optimizer (MOPSO) are also shown for performance comparisons. For the unimodal function ZDT1, the *IGD* value of the NSGWO optimization result was 0.00690, which was an order of magnitude smaller than the results of other optimization algorithms, indicating that the improved algorithm had a strong ability to jump out of the local optimum. The UF series functions were all multimodal functions, which provided many local optimal solutions, and therefore the requirements for the algorithm were relatively high. Regardless of the *IGD* value or the *GD* value, the results of the improved algorithm were better than those of the other algorithms in the UF4 function. However, the *IGD* of the improved algorithm was smaller, as was the difference between the maximum and minimum values of the *GD* value, indicating that the improved algorithm had a better stability. The above test results showed that the improved algorithm NSGWO had good stability and the ability to jump out of the local optimum.

**Table 2.** *GD* values of optimization results for benchmark functions.

**Table 3.** *IGD* values of optimization results for benchmark functions.

Figure 3 shows the Pareto fronts of optimization results of the benchmark function, in which the horizontal and vertical coordinates represent a certain target value of the benchmark function. For the ZDT1 function, the solutions of MOGWO and NSGWO were close to the real Pareto fronts. For the UF2 function, most of the MOGWO solutions were far from the real Pareto fronts, while most of the NSGWO solutions were concentrated on the real Pareto fronts. Moreover, the results indicated that the UF4 function was relatively complex. Although the optimization results of the NSGWO algorithm in Figure 3c were not close to the real Pareto fronts, they were still better than those of MOGWO. In general, the NSGWO had better accuracy and a better distribution of the obtained solutions. *Energies* **2022**, *15*, x FOR PEER REVIEW 13 of 20

**Figure 3.** Pareto fronts of multiobjective optimization algorithms for benchmark functions. (**a**) ZDT1; (**b**) UF2; (**c**) UF4. **Figure 3.** Pareto fronts of multiobjective optimization algorithms for benchmark functions. (**a**) ZDT1; (**b**) UF2; (**c**) UF4.

#### *4.2. Multiobjective Optimization of Load Dispatch Based on NSGWO 4.2. Multiobjective Optimization of Load Dispatch Based on NSGWO*

The NSGWO algorithm was employed in the multiobjective optimization of the load dispatch for coal-fired power plants while considering the economy, environmental protection, and speediness. To better show the performance of the NSGWO algorithm in load dispatch, four intelligent algorithms, including the real-number coding genetic algorithm (RCGA), particle-swarm optimizer (PSO), modified nondominated-sorting genetic algorithm (NSGA-II), and MOGWO, were employed in the following simulations. Their configurations in terms of the population size and maximum number of iterations are summarized in Table 4, and were selected through analyses of presimulations. The NSGWO algorithm was employed in the multiobjective optimization of the load dispatch for coal-fired power plants while considering the economy, environmental protection, and speediness. To better show the performance of the NSGWO algorithm in load dispatch, four intelligent algorithms, including the real-number coding genetic algorithm (RCGA), particle-swarm optimizer (PSO), modified nondominated-sorting genetic algorithm (NSGA-II), and MOGWO, were employed in the following simulations. Their configurations in terms of the population size and maximum number of iterations are summarized in Table 4, and were selected through analyses of presimulations.

**Table 4.** Configurations of NSGWO and other intelligent algorithms used in simulations of multi-

**Simulation I** 1 **Simulation II** 1 **Simulation III** 1 **NSGWO RCGA PSO NSGWO NSGA-II MOGWO NSOGWO MOGWO**

1 Simulations I-III represent the optimization of economy and environmental protection, the opti-

The NSGWO algorithm was employed for the multiobjective optimization of the load dispatch for coal-fired power plants while considering economy and environmental protection. The penalty function was used in the following simulations to deal with the equality constraints in the load-dispatch model. Firstly, a coal-fired power plant with six units was investigated; its coal consumption and pollutant-emission characteristics are listed in Table 5 [37]. The total command of these six units was 283.4 MW, and their load-changing speed was limited to 5% of the current output power. To better show the performance of the NSGWO algorithm in the load dispatch, the RCGA and PSO were introduced for a comparison with the optimization results. Then, the simulations were

4.2.1. Optimization of Economy and Environmental Protection

omy and speediness for similar-scale units, respectively.

conducted; the results are shown in Tables 6 and 7.

Population size 200 200 200 100 100 100 100 100 Number of iterations 500 500 500 300<PHONE_NUMBER> 1000

objective optimization.

**Table 4.** Configurations of NSGWO and other intelligent algorithms used in simulations of multiobjective optimization.

1 Simulations I–III represent the optimization of economy and environmental protection, the optimization of economy and speediness for various power-scale units, and the optimization of economy and speediness for similar-scale units, respectively.

#### 4.2.1. Optimization of Economy and Environmental Protection

The NSGWO algorithm was employed for the multiobjective optimization of the load dispatch for coal-fired power plants while considering economy and environmental protection. The penalty function was used in the following simulations to deal with the equality constraints in the load-dispatch model. Firstly, a coal-fired power plant with six units was investigated; its coal consumption and pollutant-emission characteristics are listed in Table 5 [37]. The total command of these six units was 283.4 MW, and their load-changing speed was limited to 5% of the current output power. To better show the performance of the NSGWO algorithm in the load dispatch, the RCGA and PSO were introduced for a comparison with the optimization results. Then, the simulations were conducted; the results are shown in Tables 6 and 7.

**Table 5.** Characteristics of dispatched units based on economy and environmental protection.

1 Variables *P* min *i* and *P* max *i* represent lower limit and upper limit of output power.

**Table 6.** Optimization results of economic load dispatch for coal-fired power plant without valvepoint effect of steam turbines.

The results showed that the performance of the NSGWO algorithm was better than that of the other algorithms in reducing power-plant coal consumption and power-plant emissions, which indicated the advantages of the NSGWO algorithm in the optimal dispatch of multiobjective economic and environmental-protection loads. Tables 6 and 7 also show that the valve-point effect of steam turbines had a great influence on the optimal dispatch results. The optimization results of Unit 1 and Unit 6 in the two models were quite different. When the valve-point effect was considered, the optimal dispatch results of these two units were 21.4 MW and 41.0 MW, respectively. However, the optimal dispatch results changed to 15.7 MW and 35.2 MW without considering the valve-point effect of the steam turbines.

**Table 7.** Optimization results of load dispatch under economy and environmental-protection strategies for coal-fired power plant considering valve-point effect of steam turbines.

In order to further illustrate the advantages of NSGWO, the Pareto fronts were determined for the NSGWO and MOGWO algorithms without considering the valve-point effect, as shown in Figure 4. The solutions of NSGWO were more widely distributed, and the solutions of MOGWO were basically dominated by the solutions of NSGWO. The comparison results show that the NSGWO algorithm had a certain application value in the multiobjective optimization of the load dispatch for coal-fired power plants. *Energies* **2022**, *15*, x FOR PEER REVIEW 15 of 20

**Figure 4.** Pareto fronts of coal-consumption cost and pollutant emissions of load-dispatch optimization under NSGWO and MOGWO algorithms. **Figure 4.** Pareto fronts of coal-consumption cost and pollutant emissions of load-dispatch optimization under NSGWO and MOGWO algorithms.

4.2.2. Optimization of Economy and Speediness for Various Power-Scale Units 4.2.2. Optimization of Economy and Speediness for Various Power-Scale Units

Considering the requirements of rapid load change of coal-fired power plants under a renewable-energy accommodation, simulations were conducted to show the performance of the NSGWO algorithm in the multiobjective optimization of economy and speediness. Here, 10 units of various power scales are selected to be dispatched in the following simulations, the characteristics of which are shown in Table 8 [38]. The total dispatch command of the 10-unit power plant was 1500 MW, and the penalty-factor values of the two-dimensional objective of economy and speediness in the NSGWO algorithm were 0.5 and 2.1, respectively. Table 9 presents the optimization results of the load dispatch under the economy and speediness strategies. To better explain the advantages of nondominated sorting, the optimization of the NSGA-II was also performed, as Considering the requirements of rapid load change of coal-fired power plants under a renewable-energy accommodation, simulations were conducted to show the performance of the NSGWO algorithm in the multiobjective optimization of economy and speediness. Here, 10 units of various power scales are selected to be dispatched in the following simulations, the characteristics of which are shown in Table 8 [38]. The total dispatch command of the 10-unit power plant was 1500 MW, and the penalty-factor values of the two-dimensional objective of economy and speediness in the NSGWO algorithm were 0.5 and 2.1, respectively. Table 9 presents the optimization results of the load dispatch under the economy and speediness strategies. To better explain the advantages of nondominated sorting, the optimization of the NSGA-II was also performed, as shown in Table 9.

**Table 8.** Characteristics of 10 dispatched units with various power scales based on economy and

11 1 1 *fP P P* ( ) 0.00510 2.2034 15 = ++ 15 7 15 60

22 2 2 *fP P P* ( ) 0.00396 1.9101 25 = ++ 20 8 20 80

33 3 3 *fP P P* ( ) 0.00393 1.8518 40 = ++ 30 10 30 100

44 4 4 *fP P P* ( ) 0.00382 1.6966 32 = ++ 35.8 15 25 120

55 5 5 *fP P P* ( ) 0.00212 1.8015 29 = ++ 50.3 8 50 150

66 6 6 *fP P P* ( ) 0.00261 1.5354 72 = ++ 75.4 6 75 280

77 7 7 *fP P P* ( ) 0.00289 1.2643 49 = ++ 124 11 120 320

88 8 8 *fP P P* ( ) 0.00148 1.2130 82 = ++ 251.85 10 125 445

99 9 9 *fP P P* ( ) 0.00127 1.1954 105 = ++ 298.75 12 250 520

10 10 10 10 *fP P P* ( ) 0.00135 1.1285 100 = ++ 298.9 9 250 550

shown in Table 9.

speediness.

1 2

2 2

3 2

4 2

5 2

6 2

7 2

8 2

9 2

10 2

**Table 8.** Characteristics of 10 dispatched units with various power scales based on economy and speediness.

**Table 9.** Optimization results of load dispatch under economy and speediness strategies for 10 units of various power scales.

1 Variables *fh* and *T* represent coal consumption per hour and adjustment time, respectively.

The results indicated that the optimal solutions of the NSGWO algorithm had the minimum adjustment times. For the minimum coal consumption of the economy strategy, the optimal solutions of the NSGWO algorithm were better than those of the other algorithms. Moreover, we found that the optimal-compromise solution of the NSGWO algorithm was able to dominate the solutions of the NSGA-II algorithm.

Figure 5 shows the Pareto fronts of the NSGWO algorithm and the MOGWO algorithm under the combined dispatch of economy and speediness. Here, a speediness index was employed to better present the optimal solutions that was defined as the square of the adjustment time. It can be seen that the optimization results of the NSGWO algorithm were more widely distributed than those of the MOGWO algorithm. The solution that dominated the optimization results of the MOGWO algorithm was found in the Pareto fronts of the NSGWO algorithm. In general, this indicated that the performance of the NSGWO algorithm in the load-dispatch optimization of economy and speediness was the best among the simulated algorithms.

was the best among the simulated algorithms.

**Figure 5.** Pareto fronts of coal consumption and speediness index of load-dispatch optimization for 10 units with various power scales under NSGWO and MOGWO algorithms. **Figure 5.** Pareto fronts of coal consumption and speediness index of load-dispatch optimization for 10 units with various power scales under NSGWO and MOGWO algorithms.

**Table 9.** Optimization results of load dispatch under economy and speediness strategies for 10

units of various power scales.

**Item Optimal-Compromise Solution Economy Strategy Speediness Strategy**

*fh*/(th −1) 1 3118.20 3120.77 3124.15 3115.02 3115.22 3123.25 3130.11 3131.29 3126.27 *T*/min 1 4.91 4.98 4.49 7.48 7.12 5.13 3.52 3.59 3.77 *P1*/MW 15.64 20.96 24.48 15.14 15.00 23.50 34.21 35.94 25.66 *P2*/MW 37.58 37.88 45.04 26.00 27.33 45.24 44.566 45.50 46.30 *P3*/MW 45.16 48.03 48.21 36.20 34.21 45.48 59.17 61.00 52.29 *P4*/MW 65.82 69.80 76.64 58.34 55.86 76.33 80.19 70.07 77.91 *P5*/MW 79.48 76.81 76.43 76.74 82.85 75.24 77.36 78.82 78.27 *P6*/MW 101.45 94.49 97.10 109.45 109.63 95.54 94.91 96.26 97.04 *P7*/MW 159.34 156.95 152.18 146.84 149.06 152.97 157.04 155.41 156.39 *P8*/MW 298.22 294.80 296.78 302.62 306.39 301.44 284.90 286.40 289.48 *P9*/MW 354.25 358.46 345.97 364.46 356.67 339.19 339.36 341.79 344.03 *P10*/MW 343.06 341.81 337.18 364.15 363.01 345.06 328.23 328.81 332.64

**NSGWO NSGA-II MOGWO NSGWO NSGA-II MOGWO NSGWO NSGA-II MOGWO**

1 Variables *fh* and *T* represent coal consumption per hour and adjustment time, respectively.

algorithm was able to dominate the solutions of the NSGA-II algorithm.

The results indicated that the optimal solutions of the NSGWO algorithm had the minimum adjustment times. For the minimum coal consumption of the economy strategy, the optimal solutions of the NSGWO algorithm were better than those of the other algorithms. Moreover, we found that the optimal-compromise solution of the NSGWO

Figure 5 shows the Pareto fronts of the NSGWO algorithm and the MOGWO algorithm under the combined dispatch of economy and speediness. Here, a speediness index was employed to better present the optimal solutions that was defined as the square of the adjustment time. It can be seen that the optimization results of the NSGWO algorithm were more widely distributed than those of the MOGWO algorithm. The solution that dominated the optimization results of the MOGWO algorithm was found in the Pareto fronts of the NSGWO algorithm. In general, this indicated that the performance of the NSGWO algorithm in the load-dispatch optimization of economy and speediness

#### 4.2.3. Optimization of Economy and Speediness for Similar-Scale Units

Moreover, a real coal-fired power plant with eight 600 MW units and two 660 MW units was investigated under a renewable-energy accommodation, which may further explain the performance of the NSGWO algorithm in the load dispatch. These 10 units had similar power scales, so a higher optimization accuracy was needed in the algorithm. The coal-consumption characteristics, as well as other initial conditions and constraints, are shown in Table 10, in which the unit for coal consumption is grams of coal per kWh. In addition, the penalty-factor values of the two-dimensional targets of economy and speediness in the algorithms were 0.5 and 1.1, respectively, and the total power command of this 10-unit power plant was 4000 MW.

**Table 10.** Characteristics of 10 dispatched units based on economy and speediness with 10 similarscale units.

Table 11 shows the optimization results for the real coal-fired power plant with 10 similar-scale units. It indicates that the minimum coal consumption of NSGWO was 328 g/(kWh), which was better than that of MOGWO. The minimum adjustment time of NSGWO was 1.36 min, which was also shorter than that of MOGWO. Furthermore, the solutions of the MOGWO algorithm for these two minima were dominated by the NSGWO algorithm, which meant that the distribution of the solutions of the NSGWO algorithm was better than that of the MOGWO algorithm.

**Table 11.** Optimization results of load dispatch under economy and speediness strategies for a real coal-fired power plant with 10 similar-scale units.

1 Variables *fw* and *T* represent coal consumption per kWh and adjustment time, respectively.

Figure 6 illustrates the Pareto fronts of coal consumption and adjustment time of the former simulations for NSGWO and MOGWO. It shows that the optimization results of NSGWO were basically at the lower left of MOGWO, indicating that the optimization effect of NSGWO was better. Most of the solutions of the MOGWO algorithm are concentrated in the upper left, while the solutions of the NSGWO algorithm are more uniformly distributed. This indicates that the solutions of the NSGWO algorithm had a better distribution in the combined load dispatch while considering the requirements of both economy and speediness for coal-fired power plants. Figure 6 illustrates the Pareto fronts of coal consumption and adjustment time of the former simulations for NSGWO and MOGWO. It shows that the optimization results of NSGWO were basically at the lower left of MOGWO, indicating that the optimization effect of NSGWO was better. Most of the solutions of the MOGWO algorithm are concentrated in the upper left, while the solutions of the NSGWO algorithm are more uniformly distributed. This indicates that the solutions of the NSGWO algorithm had a better distribution in the combined load dispatch while considering the requirements of both economy and speediness for coal-fired power plants.

**Figure 6.** Pareto fronts of coal consumption and adjustment time of load-dispatch optimization for coal-fired power plant with 10 similar-scale units under NSGWO and MOGWO algorithms. **Figure 6.** Pareto fronts of coal consumption and adjustment time of load-dispatch optimization for coal-fired power plant with 10 similar-scale units under NSGWO and MOGWO algorithms.

The results of this study showed that the proposed nondominated-sorting GWO algorithm achieved high accuracy and had a good distribution of solutions, which may extend its applications in multiobjective optimization. We caution that the presented analyses for the NSGWO and other intelligent algorithms were somewhat simplified, and further investigations are needed. The optimization procedure proposed in this study may help researchers conduct other investigations into the optimization of multiobjective problems and the operation of power-generation systems. The results of this study showed that the proposed nondominated-sorting GWO algorithm achieved high accuracy and had a good distribution of solutions, which may extend its applications in multiobjective optimization. We caution that the presented analyses for the NSGWO and other intelligent algorithms were somewhat simplified, and further investigations are needed. The optimization procedure proposed in this study may help researchers conduct other investigations into the optimization of multiobjective problems and the operation of power-generation systems.

dation were conducted based on the NSGWO algorithm. On the basis of the performed

(1) The proposed NSGWO algorithm employed efficient nondominated sorting, a reference-point selection strategy, and a simulated binary crossover operator in order to avoid falling into the local optimum, enhance the solution diversity, and shorten the convergence time. The optimization results of benchmark functions indicated that the NSGWO algorithm had a better accuracy and a better distribution than the

(2) Regarding the load dispatch of economy and environmental-protection strategies, the performance of the NSGWO algorithm was better than those of the other simulated algorithms in reducing power-plant coal consumption and emissions. In addition, the widely distributed solutions of the NSGWO algorithm could dominate those of the MOGWO algorithm, which elucidated the application value of the NSGWO algorithm in the multiobjective load dispatch of coal-fired power plants.

analyses, the conclusions of this paper can be summarized as follows:

**5. Conclusions**

MOGWO algorithm.

#### **5. Conclusions**

Multiobjective optimizations of load dispatch under a renewable energy accommodation were conducted based on the NSGWO algorithm. On the basis of the performed analyses, the conclusions of this paper can be summarized as follows:

**Author Contributions:** Conceptualization, Y.C. and F.S.; methodology, Y.C. and T.L.; software, T.L. and T.H.; validation, T.H.; formal analysis, Y.W. and M.L.; investigation, Y.C. and T.L.; writing original draft preparation, Y.C.; writing—review and editing, T.H. and F.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Natural Science Foundation of Jiangsu Province (grant number BK20210240).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available upon request from the corresponding author.

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

#### **References**

## *Article* **Reserve-Constrained Unit Commitment Considering Adjustable-Speed Pumped-Storage Hydropower and Its Economic Effect in Korean Power System**

**Woo-Jung Kim 1 , Yu-Seok Lee 1 , Yeong-Han Chun 1,\* and Hae-Seong Jeong 2**

**Abstract:** The Korean government has declared the goal of net-zero-carbon emissions with a focus on renewable energy expansion. However, a high proportion of baseload generators and an increasing proportion of variable renewable energy (VRE) may cause problems in the power system operation owing to the low cycling capability of baseload generators and variability of VRE. To maintain system reliability, the government is planning to construct pumped-storage hydropower (PSH) plants, which can provide flexibility to the system. This study evaluated the operating cost savings obtained by different types of PSH: the adjustable-speed PSH (AS-PSH) and fixed-speed PSH (FS-PSH), based on the duck-curve phenomenon and the increase in spinning reserve requirement. In this study, the reserve-constrained unit commitment was formulated using a mixed-integer-programming considering the operational characteristics of AS-PSH and conventional generators. To consider the duck-shaped net-load environment, the projected VRE output data were calculated through physical models of wind turbines and photovoltaic modules. The operating costs for the non-PSH, FS-PSH, and AS-PSH construction scenarios were KRW 43,129.38, 40,038.44, and 34,030.46, respectively. The main factor that derived this difference was determined to be the primary reserve of AS-PSH's pumping mode.

**Keywords:** adjustable-speed pumped-storage hydropower; fixed-speed pumped-storage hydropower; variable renewable energy; duck-curve; reserve-constrained unit commitment; spinning reserve requirement; nuclear power generators

## **1. Introduction**

## *1.1. Research Motivation*

The Korean government has declared the goal of achieving net-zero carbon emissions by 2050. Therefore, Korea is planning to significantly increase the proportion of renewable energy resources. In the generation expansion plan, which determines the Korean investment in power generation facilities, 2030 was set as the target year to ensure that 20% of the annual power is generated by renewable energy resources [1]. Among the renewable energy resources, 75% is expected to be generated from wind and solar energy, which are variable renewable energy (VRE) resources. However, VRE resources may cause difficulties in balancing the supply and demand within a power system owing to their output variability and uncertainty [2]. Specifically, solar power generation makes the system net-load have the characteristics of a duck-curve, which is a concept first introduced by the California Independent System Operator [3]. It has a duck-shaped net-load pattern, which is carried by conventional generators. The VRE's long-term and short-term variability require increased cycling capabilities of the generator and operating reserves, respectively. Therefore, it is difficult for the system operator to maintain system reliability while maintaining the operating proportion of baseload generators during the daytime.

**Citation:** Kim, W.-J.; Lee, Y.-S.; Chun, Y.-H.; Jeong, H.-S. Reserve-Constrained Unit Commitment Considering Adjustable-Speed Pumped-Storage Hydropower and Its Economic Effect in Korean Power System. *Energies* **2022**, *15*, 2386. https://doi.org/ 10.3390/en15072386

Academic Editor: Massimiliano Renzi

Received: 3 January 2022 Accepted: 18 March 2022 Published: 24 March 2022

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

**Copyright:** © 2022 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/).

In 2020, the first output reduction of nuclear power generators was conducted by the system operator to maintain reliability in the Korean power system. Shinkori nuclear power generators (#3 and #4) have a total capacity of 2800 MW. The power system operator was scheduled to reduce the total output power to 600 MW for 8 h, and the scheduled output power was maintained for 13 h. The system operator predicted that the system demand was lower than usual owing to long special holidays and analyzed that it would be difficult to meet the system frequency maintenance range if a nuclear power generator with the maximum output (1400 MW) is tripped. There are similar issues with implementing net-zero carbon emissions. Although the coal-fired and nuclear power generators are expected to gradually phase out, the proportion of these baseload generators is still expected to be high in 2030. It is especially difficult to expect frequent output adjustments or startup/shutdown operations with nuclear power generators. In addition, existing nuclear power generators in Korea do not provide frequency response capabilities. Accordingly, situations can occur where several nuclear power generators would need to be operated at a minimum generation level or even excluded from scheduling. In those cases, the operating cost of the power system can increase significantly, and the existing baseload generators may even become stranded assets.

Grid-scale energy storage systems, such as battery energy storage systems (BESS) and pumped-storage hydropower (PSH), can provide solutions to operating costs by mitigating the effect of the duck-curve phenomenon of the net-load and increasing the utilization of nuclear power generators. Although BESSs are capable of providing load shifting and fast frequency responses, PSH is technically proven and more desirable in that it has a large energy capacity and the capability of providing operating reserves. The Korean government is planning to build PSH plants [1]. However, there are three types of PSH plants: fixed-speed (FS), adjustable-speed (AS), and ternary (T) types. As a quickstart resource, a conventional FS-PSH can be rapidly deployed to the system in case of a generator's forced outage or insufficient reserve [4]. AS-PSHs and T-PSHs can provide the same operational advantages as FS-PSHs. In addition, AS-PSH can control its rotor speed using power electronic devices, which enables a rapid power adjustment in the generating and pumping modes and is therefore expected to further improve the flexibility of the power system. Given that the pump and turbine of T-PSH can rotate simultaneously on the same shaft in the same direction, it has the advantage of adjusting the pumping power in the hydraulic short-circuit mode and requires less time to transition between the generating and pumping modes [5]. Therefore, we needed to decide which type will be of benefit to our system when they are constructed, which was the motivation of our research. In this paper, we focused on the economic effect obtained by AS-PSH using the unit commitment (UC) optimization approach in a future Korean power system. In the case of T-PSH, there are three modes: generating mode, pumping mode, and short-circuit mode, which each have different operational characteristics. Therefore, future research is needed to model the operational characteristics of T-PSH in UC formulation.

#### *1.2. Related Research*

Typically, PSH charges at night when the system demand is low and discharges during peak-demand hours. However, owing to the large-scale integration of VREs, the economic benefits that can be obtained from the PSH have been studied from new perspectives [6–12]. In [6], the impact of reducing the operating cost achieved by the input of FS-PSH in a power system with high wind energy was analyzed, and FS-PSH was found to reduce the scheduling cost by reducing the startup and shutdown costs of thermal power generators. In [7], the effect of reducing operating costs through FS-PSH in a wind energy environment was analyzed through a genetic algorithm. The paper derived the reduction in energy purchase cost and environmental cost of the power system. In [8], it was found that as the initial upper reservoir volume of PSH increases, the overall operating cost consisting of grid power purchase cost, demand response program purchase cost, and solar power sale income decreases. In [9], the effect of reduction in the wind power curtailment and the

total operating cost consisting of energy cost, ancillary service cost, and startup cost were derived using FS-PSH as a reserve resource. In [10], the energy cost in a standalone hybrid wind/solar/biomass system was lower in the scenario when PSH was the input than when battery storage was the input. In [11], the authors evaluated the potential contribution of PSH to the Greek power system's future scenario. Through UC simulation, it was deduced that the PSH contributes to energy balance by absorbing the over-generation of wind and solar energy. Ref. [12] proposes a stochastic scheduling approach considering price-based demand response and PSH in a microgrid. In [12], the operating cost was reduced by PSH's discharging at the peak time and charging at the off-peak time. However, given the proposed methods in [6–12], it is difficult to verify the economic contribution obtained from the flexibility of AS-PSH.

The reduction in operating costs due to the flexibility of AS-PSH was studied in [13–15]. In [13], the dynamic dispatch optimization problem was solved with the reserve and network constrained in a wind energy environment, and an economic analysis was performed on FS-PSH, T-PSH, and AS-PSH. However, in [13], it was difficult to verify probable changes in operating conditions within the duck-curve phenomenon, given that the influence on solar energy was not considered. In [14], the economic impact of AS-PSH was analyzed under the conditions of minimizing fuel, cycling, and startup cost. It was modeled as having an operating range in the pumping mode and was derived to reduce the operating cost. In [15], the cost benefits obtained with AS-PSH from the viewpoint of daily operation were compared with those of FS-PSH. However, the study focused on the AS-PSH's flexible operating range in pumping. Therefore, it is still necessary to focus on the reserve capability of AS-PSH that can respond to VRE's short-term variations.

In addition to PSH-related studies, UC optimization research is ongoing from various perspectives. Authors in [16–18] developed methods to efficiently find the UC solution at a high time resolution to cope with the variability and intermittent nature of renewable energy. Ref. [16] proposes a thermal unit commitment algorithm through dynamic programming that combines a priority list method to perform UC and ED at a time resolution of less than 60 min. In [17], the authors pointed out that the branch and cut strategy may exceed the allowable value in the sub-hourly UC problem due to large numbers of virtual transactions. In [17], the authors reduced the calculation time by including the concept of ordinal-optimization in surrogate absolute-value Lagrangian relaxation [19]. Ref. [18] proposes an approach to reduce calculation time by implementing soft constraints of reserve and transmission capacity within the surrogate absolute-value Lagrangian relaxation framework. References [20,21] analyzed the impact of greenhouse gas emissions on the UC problem. In [20], the authors formulated the environmental UC problem considering greenhouse gas emission constraints and carbon taxes in the existing UC approach. In [21], the authors adopted the piecewise linearized greenhouse gas function of the generator's loading level to the UC based on the MIP method. The greenhouse gas function was applied in the form of cost by applying a weighting factor.

#### *1.3. Contribution*

Previous studies on PSH verified its effects on the operating cost savings of power systems with integrated VRE and PSH. However, the majority of these studies focused on the traditional FS-PSH's contribution to power systems in terms of energy saving. The other studies [13–15] on AS-PSH did not consider its capability to supply primary and secondary reserves, particularly in the pumping mode. In a system with a considerable amount of solar energy resources, system operators need large operating reserves, and most generation sources need to be off-lined owing to an over-generation; on the other hand, most PSHs are required to operate in pumping mode. Therefore, to analyze the economic contribution of AS-PSH in a power system with a high proportion of renewable energy resources, it is necessary to model the AS-PSH to provide various reserves, including the primary reserve. Based on the above discussion, the purpose of this study is to verify the

economic contribution of AS-PSH's reserve capability compared with that of FS-PSH in the Korean power system.

To achieve the objective of this study, a reserve-constrained UC (RCUC) was formulated using mixed-integer-programming (MIP) to consider the cycling characteristics of conventional thermal generators in the duck-curve environment and reflect the Korean operating reserve services. The primary, secondary, and regulation reserves defined by the Korean power market were modeled in the RCUC formulation. This formulation can be widely used in other power systems with a similar scheme of system frequency restoration by primary and secondary reserves. The time-series data of the VRE were calculated through physical models, including the output model of wind turbines and photovoltaic (PV) modules. The short-term variation of the VRE was extracted from the time-series data. The increment in the spinning reserves was considered under the n-sigma criterion, assuming that the distribution of the VRE's short-term output variation follows a Gaussian distribution. With the operation range set in the generating and pumping modes, the AS-PSH was modeled to provide primary, secondary, and regulation reserves. Simulations were performed on the 2030 generation mix in Korea. Scenarios were set depending on the types of new PSH plants, namely, non-PSH, FS-PSH, and AS-PSH plants. The simulation was implemented using the optimization solver Xpress-MP.

The main contributions in this paper are listed below:

#### **2. Variable Renewable Energy in Future Korean Power System**

Before analyzing the influence of the new PSH on the future power system in Korea, a methodology to generate the time-series data of the 2030 VRE is introduced. Based on the data, the increment in the spinning reserve due to the short-term variation of the VRE was estimated.

#### *2.1. Modeling of Variable Renewable Energy's Output Power*

The 2030 target capacity of the VRE is about 54.2 GW, and wind and solar energy account for 17.7 GW and 36.5 GW, respectively [1]. Based on the capacity of the VRE to be installed, the 2030 output power can be calculated using historic regional meteorological data. From the Korea Meteorological Administration (KMA), the regional wind speed, solar irradiance, and temperature data measured every minute were used [22].

#### 2.1.1. Wind Energy Output Power Model

To calculate the total output power of wind energy in 2030, it is necessary to determine the location and capacity of each wind farm. To solve the problem, it was assumed that the total capacity of wind energy in which the location is not determined would be distributed to areas with high energy potentials. Therefore, the capacity was distributed proportionally to the locations with high average annual wind speed. Finally, a total of 94 locations of the 2030 wind farms and their capacity were determined.

In a wind farm, each wind turbine is installed at a distance. Therefore, input data conversion was performed without applying the same wind speed to all wind turbines in the same wind farm. Each wind farm consists of N turbine units, and the wind speed after (*n* − 1) minutes is applied to calculate the output power of the nth wind turbine at time *t*.

$$P\_{\rm WF}(t) = \sum\_{n=1}^{N} P\_{\rm WT,n}(t) = \sum\_{n=1}^{N} P\_{\rm WT,n}(V\_h(t - n + 1))\tag{1}$$

To calculate the output power of the wind turbine, Equations (2) and (3) were used [23]. For the convenience of calculation, each wind turbine group was modeled with one unit. The data in Table 1 were used for the specifications of onshore and offshore wind turbines [24].

$$P\_{WT,n}(t) = \begin{cases} 0 & \text{for } V\_h(t-n+1) < V\_{\text{in}} \text{ and } V\_h(t-n+1) > V\_{\text{out}}\\ P\_{\text{curv}} \text{for}(V\_h(t-n+1)) \cdot \frac{\text{C}\_{\text{NF}}}{\text{P}\_r \cdot \text{N}} & \text{for } V\_{\text{in}} \le V\_h(t-n+1) \le V\_{\text{r}}\\ \frac{\text{C}\_{\text{NF}}}{\text{N}} & \text{for } V\_r < V\_{\text{ll}}(t-n+1) \le V\_{\text{out}} \end{cases} \tag{2}$$

$$P\_{\rm current}(V\_{\hbar}(t)) = a\_1 V\_{\hbar}(t)^6 + a\_2 V\_{\hbar}(t)^5 + a\_3 V\_{\hbar}(t)^4 + a\_4 V\_{\hbar}(t)^3 + a\_5 V\_{\hbar}(t)^2 + a\_6 V\_{\hbar}(t) + a\_7 \tag{3}$$

**Table 1.** Parameters of wind turbines (onshore and offshore).

Given that wind speed increases with an increase in altitude, the wind speed used as input data for Equations (2) and (3) should be measured at the tower height. However, the wind speed data of the KMA are measured at an altitude of 10 m. Therefore, using Equation (4) based on the wind power law, we estimated wind speed data at the tower height [25].

$$V\_h = V\_m \left(\frac{Z\_h}{Z\_m}\right)^a \tag{4}$$

2.1.2. Solar Energy Output Power Model

By 2030, 36.5 GW of PV plants is expected to be installed. As in the case of wind power, it is necessary to determine the capacity and the location of each PV plant. However, most PV plants are expected to be installed at a small scale and dispersed over the nation. Therefore, it is necessary to estimate insolation data for areas where KMA does not observe the insolation. To address this problem, we expanded the number of observation stations using a two-dimensional interpolation method [26]. The capacity of each PV plant was assigned in proportion to historical accumulated solar insolation.

The solar insolation data were converted to the irradiance data. The PV output power, which is determined by the irradiance and temperature on the PV module surface, can be expressed as Equation (5) [27].

$$P\_{PV,M} = P\_{PV,M}^\* \cdot \frac{G}{G^\*} [1 - \gamma \cdot (T\_{\mathbb{C}} - 25)] \tag{5}$$

Given that each PV plant consists of PV arrays of several PV modules, the output power of each PV plant can be calculated using Equation (6).

$$P\_{PV} = P\_{PV,M} \cdot \mathbb{C}\_{PV} / P\_{PV,M}^\* \tag{6}$$

Given that the temperature on the PV module surface used in Equation (5) is difficult to measure, it can be approximated using Equation (7) based on the irradiance and the air temperature [28]. to measure, it can be approximated using Equation (7) based on the irradiance and the air temperature [28]. , ∙/∗(6)Given that the on the PV module surface in Equation is to measure, it can approximated using Equation on the irradiance and the air temperature [28].

Given that each PV plant consists of PV arrays of several PV modules, the output

∙ ∗ [1 − ∙ − 25)] (5)Given that each plant consists of PV arrays of several PV modules, the output

∗

= , ∙ /,

Given that the temperature on the PV module surface used in Equation (5) is difficult

$$T\_{\mathbb{C}} = T\_A + \frac{T\_{\text{NOCT}} - 20}{800}G\tag{7}$$

[1 − ∙ ( − 25)] (5)

(6)

#### *2.2. Characteristics of Variable Renewable Energy in Korea 2.2. Characteristics of Variable Renewable Energy in Korea 2.2. Characteristics of Variable Renewable in Korea*

*Energies* **2022**, *15*, x FOR PEER REVIEW 6 of 24

power of each PV plant can be calculated using Equation (6).

power of each plant can calculated using Equation (6).

, = , ∗ ∙ ∗

= , ∗

REVIEW of

In this subsection, the output power characteristic of the 2030 VRE is described. Figure 1 depicts the calculated time-series output power of the VRE and its duration curve for the year 2030. In this subsection, the output power characteristic of the 2030 VRE is described. Figure 1 depicts the calculated time-series output power of the VRE and its duration curve for the year 2030. In this the output of the 2030 VRE ure 1 depicts the time-series output of the VRE and its curve for the year 2030.

**Figure 1.** Output and duration curve of VRE for 2030. **Figure 1.** Output and duration curve of VRE for 2030. **Figure** Output and curve of

Figure 2 presents the average daily wind and solar power generation of each month. The average daily wind power generation is highest in winter and lowest in summer. The average daily solar power generation is highest in spring and lowest in winter. Based on these results, the average daily VRE power generation is highest in spring and lowest in autumn. Figure 2 presents the average daily wind and solar power generation of each month. The average daily wind power generation is highest in winter and lowest in summer. The average daily solar power generation is highest in spring and lowest in winter. Based on these results, the average daily VRE power generation is highest in spring and lowest in autumn. Figure 2 presents average each month. The average daily power generation highest in and lowest in The average generation is highest in and lowest in Based on these results, average VRE generation is highest spring lowest in autumn.

**Figure 2.** Average daily VRE power generation of each month for 2030. **Figure** daily VRE generation of month for 2030. **Figure 2.** Average daily VRE power generation of each month for 2030.

Figure 3 illustrates the hourly average output power for wind and solar energy throughout the year 2030. Wind energy demonstrates the highest output power between 14:00 and 15:00, and the highest PV output power and both wind and solar output power are between 12:00 and 13:00.

are between 12:00 and 13:00.

**Figure 3***.* Hourly average output (MW) of wind and PV power projected for 2030. **Figure 3.** Hourly average output (MW) of wind and PV power projected for 2030.

*2.3. Influence of Variable Renewable Energy on Reserve Requirements of the Korean Power System*

Figure 3 illustrates the hourly average output power for wind and solar energy throughout the year 2030. Wind energy demonstrates the highest output power between 14:00 and 15:00, and the highest PV output power and both wind and solar output power

*2.3. Influence of Variable Renewable Energy on Reserve Requirements of the Korean Power System* Traditionally, most power system operators have used the N−1 contingency method, which secures a reserve equal to one or more largest units [29]. The system operator of the Korean power system uses the traditional reserve procurement approach to cover the larg-Traditionally, most power system operators have used the N−1 contingency method, which secures a reserve equal to one or more largest units [29]. The system operator of the Korean power system uses the traditional reserve procurement approach to cover the largest generator outage (1400 MW). Table 2 outlines the operating reserve requirements for the Korean power market [30].

est generator outage (1400 MW). Table 2 outlines the operating reserve requirements for

the Korean power market [30]. **Table 2.** Operating reserve requirements for the Korean power system.

Quick-start reserve 2000 MW Standstill Manual Load forecast error, etc.

The primary reserve prevents the system frequency from dropping immediately after an instantaneous event, such as a generator's forced outage by the synchronized generator's governor response, and the secondary reserve is used to restore the system frequency to the nominal frequency (60 Hz). The required primary reserve in Table 2 was calculated based on the value at the time of a quasi-steady state of the system frequency after the outage of the largest unit (1400 MW) in Korea. In this process, an unused capacity of 1400 MW is required as the secondary reserve to restore the system frequency to the nominal frequency (60 Hz). Immediately after restoration to the nominal frequency by the secondary reserve, the primary reserve is automatically resecured. The regulation reserve is estimated as 700 MW to respond to non-instantaneous events, such as short-term load fluctuations. The tertiary reserve is secured by resources in the standstill state to cover the The primary reserve prevents the system frequency from dropping immediately after an instantaneous event, such as a generator's forced outage by the synchronized generator's governor response, and the secondary reserve is used to restore the system frequency to the nominal frequency (60 Hz). The required primary reserve in Table 2 was calculated based on the value at the time of a quasi-steady state of the system frequency after the outage of the largest unit (1400 MW) in Korea. In this process, an unused capacity of 1400 MW is required as the secondary reserve to restore the system frequency to the nominal frequency (60 Hz). Immediately after restoration to the nominal frequency by the secondary reserve, the primary reserve is automatically resecured. The regulation reserve is estimated as 700 MW to respond to non-instantaneous events, such as short-term load fluctuations. The tertiary reserve is secured by resources in the standstill state to cover the case of reserve shortages, and the quick-start reserve responds to load forecasting errors.