Optimal Operation of a Photovoltaic Integrated Captive Cogeneration Plant with a Utility Grid Using Optimization and Machine Learning Prediction Methods

: The World Energy Council, in its 2019 World Energy Scenarios Report, advised policy-makers to identify innovative opportunities for the integration of renewable energy resources into existing electrical power systems to achieve a fast and affordable solution. However, large-scale industries with cogeneration units are facing problems in handling the higher penetration levels of intermittent renewable energies. This paper addresses large-size photovoltaic power integration problems and their optimal operation. This work considers the case of a chemical industry having both cogeneration power and solar photovoltaics. Here, a modiﬁed ﬁreﬂy algorithm and a hybrid power resource optimization solver are proposed. The results of the proposed method are compared with other benchmark techniques, to conﬁrm its advantages. The proposed techniques can be used in industries having cogeneration power plants with photovoltaics for better optimization and to meet the guidelines speciﬁed in IEEE 1547. The voltage ramp index is proposed to determine the voltage ramp up and down with intermittent solar irradiance. Additionally, a machine learning technique is used to predict the cogeneration plant efﬁciency at different loads and the solar irradiance under varying weather conditions. Finally, this paper proposes the effectiveness of the modiﬁed heuristic technique and certain guidelines, including solvers for industrial use. The linear regression–linear model results of the RMSE of 0.2313, R -squared of 0.98, MSE of 0.05351, MAE of 0.2018 and p -values of 5.2336 × 10 for the and that the regression model is ﬁt for ﬁeld CGPP is suitable controllers for its optimum operations, up to the lowest acceptable of operations operate efﬁciently and are not advisable.


Introduction
The 2019 World Energy Scenarios Report has predicted an increase in the share of renewable power from the present 26 % of total power generation in 2020 to 33% in 2040 [1]. The radical shift in the electrical power industry is due to the three 'D's of decarburization (with renewable energy resources), decentralization (with distributed generation) and digitalization (artificial intelligence) [2]. First, decarburization is demanding for harnessing renewable energy resources (RERs) such as solar, wind, geothermal, tidal, etc. The traditional vertically designed power systems are being replaced with the latest distribution generations. Second, this type of shift from centralized generation to decentralized generation or distributed generation (DG) has made the existing distribution networks more complicated due to two-way power flows. Various benefits and challenges, including their classification, are mentioned in [3]. The classification is mainly based on the type of technology such as renewables (PV, wind, etc.) and non-renewables (gas turbines, reciprocating engines, etc.). Additionally, there are other classifications based on their size, and some of them only inject active power. Hence, there is a need for proper integration of DGs into existing distribution networks. This also demands standard microgrid structures and controllers [4]. Third, the evolution of digital technologies has made complex operations somewhat simpler. Digital tools are helping in better management and optimization of various power resources, especially RERs. The artificial intelligence (AI) techniques such as machine learning (ML) are emerging in the power sector for better RE prediction, which is mostly intermittent. Industries should be able to predict hour-or a day-ahead power with varying weather conditions to rely on RERs for their optimization and better control. For this purpose, a more useful solution with machine learning is also proposed.
The optimal operation of industries with cogeneration facilities is becoming difficult due to the integration of solar photovoltaic (SPV) power. The major integration problems are maintaining the stability, voltage profile and efficiency of cogeneration power plant (CGPP) generators with intermittent PV power. Additionally, maintaining the synchronism and applicable codes of utility grids is a major problem of concern. Various heuristic techniques exist for the optimal utilization of highly stochastic RERs with a special mention of SPV. However, little work on optimum power allocation problems in industries with cogeneration power plants has been performed. Hence, there is a need for a modified heuristic to cater to the needs of industries having CGPPs with SPP integration. With the above-cited practical insights and evolution of AI techniques in electrical power systems, the authors have collected data from a large chemical process industry having both CGPP and SPP and conducted the study on the impacts of integration and the operational strategies of the industry. With these motivations from this problem, a heuristic technique for the optimal operation of a cogeneration plant is developed.

Literature Survey
A detailed literature survey was conducted to determine the RE integration problems, the various techniques adopted and the gaps therein. The installation of RERs in the form of DG in the existing power distribution network may be attributed to the reduced line losses and enhanced power transfer capacity. Conventional and metaheuristic algorithms are being used to solve short-term power scheduling, such as hydrothermal power scheduling, but sometimes they may lead to premature convergence [5]. Additionally, IEEE 1547-2018 [6] with its amendment-1, noted as 1547a-2020 [7], has given certain guidelines on RE integration but not on cogeneration integration problems. According to the 2/3 rule [8] in a radial feeder, 2/3 of the capacity of the incoming distributed capacity should be installed at 2/3 of the length of the distribution system. However, this has certain operational impacts such as regularity requirements, voltage profiles, and operational modes. However, their optimum utilization along with existing utility grids and cogeneration power systems can be ensured using either analytic or heuristic techniques. Optimization techniques are mainly classified into analytical (such as the Newton-Raphson method and linear programming) and heuristic methods (such as particle swarm optimization). The authors in [9] mentioned that energy management should consider operational constraints and suggested the optimal utilization of PV power and other resources. The authors suggested caution due to the premature failure of electrical equipment due to rapid changes in RE power and suggested short-term forecasting of PV power and stochastic scheduling of various power resources. The impact of RESs as cogeneration power plants on power system networks is discussed with an example of a 6.3 MW cogeneration power plant in [10]. This work was performed using a small cogeneration plant and suggested managing the reactive power requirements. The authors in [11] reviewed the integration of PV into electricity grids and suggested grid forecasting to improve the grid's health in real-time and timely intervention with intelligent technologies to correct the problems. An expert group of the International Energy Agency (IEA) on recommended practices for wind and PV integration [12] provided information on how to perform integration studies with main points such as load scheduling and dispatch, reactive power and transients for future studies. Furthermore, they mentioned that RE integration methods would continue to evolve with real-time studies and suggested the development of flexible tools for use in different time periods.
The authors in [13] reviewed the RE policies in India and the barriers to PV generation such as lack of awareness and an absence of adequate energy management systems. To enable reactive power capability in PV inverters, the authors in [14] suggested installing Energies 2021, 14, 4935 3 of 28 suitable capacitors when the power factor of the SPP is assumed to be unity. Additionally, the reactive power capability of PV inverters can be actively controlled without incurring much burden on the inverters, subsequently increasing the hosting capacity of renewable energy in microgrids [15]. The authors in [16] performed a case study on the reactive power management of a PV integrated cogeneration power plant to determine the effects of higher PV integrations and protection issues [17] and suggested certain guidelines for RE integration. The authors in [18] have studied the challenges with high RE penetration of Croatia and confirmed that even a 30% RE penetration itself will cause voltage and frequency problems. They have suggested to adopt novel methodologies to mitigate the higher RE penetration effects. The authors of [19] examined grid operators' experience in handling high penetration of RERs of up to 100% and suggested some new technologies for such higher penetrations. The effects of PV penetration on a distribution system with 100 feeders were statistically studied in [20], and it was concluded that 20% PV penetration requires strengthening the existing system to avoid voltage problems. The work also recommended conducting field studies as further work in this area. In [21], a prototype microgrid with PV, wind and a battery storage system seeking to achieve proper energy management was developed. In this experiment, the authors conducted tests by varying irradiance and wind speed levels without knowing the actual weather conditions. The authors in [22] also reviewed the effects of PV integration on both active and reactive powers. They have also discussed the curtailment of active power during over-voltages and reactive power management during low voltage conditions and suggested PV forecasting for voltage stability management. In [23] the authors reviewed various machine learning techniques such as regression analysis and metaheuristic techniques for PV power forecasting. It was mentioned that this work can be used for future improvements. The authors in [24] discussed the higher RER penetration effects, introduced a hybrid irradiance forecasting technique and mentioned further improvements in the future using past field data.
The authors of [25] reviewed various RER integration models, optimization techniques, etc., and introduced the monarch butterfly algorithm. After discussing various algorithms such as the GA and PSO, the authors stressed the need to address the modeling of the uncertainty of RERs, the load demand, etc. The authors in [26] discussed the energy transition to 2050 and projected the share of RERs at 63% compared to 15% in 2015. They have indicated the need for new technological innovations and new business models. Further research and the findings are mentioned in the respective sections. The literature review indicates that there are many gaps in the existing system, and there is a need to bridge the gaps with new technologies, optimization and integration methods. Hence, there is a need to have a more suitable optimization method for hybrid resources, especially industries with cogeneration power plants integrated with large solar PV plants. It is also required to validate the power prediction and energy optimization techniques, including the technical validation of optimization results with any well-established power system analysis software.
The contributions of this work are given below: • All the components of the power distribution network were modeled with LinDistFlow branch equations.

•
The field data and machine learning were used to find the relation function between the percentage load to forecast the CGPP efficiency for field operator use, and field data were used for validation. • Field and online available weather data and linear regression analysis were used to obtain the relation function between site weather conditions and irradiance. The same has been validated with field data, and the method can be used to predict site-specific time PV generation for varying weather conditions to achieve better load management. • A hybrid power resources optimization solver for industrial use, which can avoid the expert requirement of GAs, was developed and demonstrated. • Firefly algorithm (FA), a metaheuristic algorithm that generally has better performance than other methods, was implemented in the operation of industries with cogeneration plants having PV power and connected with utility grids. • A modified firefly algorithm (MFA) for the above-mentioned industries was developed and demonstrated, and it achieved further improved performance when compared with the base case of the FA.

•
The results obtained were technically validated using the solver with E-Tap software.

•
The integration guidelines given in the standards such as IEEE 1547-2018 with its amendment 1as 1547a-2020 were fully followed.
The remainder of this work is organized in the following way. Starting with the literature survey in Section 2, brief details of the industry under study are mentioned in Section 3. Sections 4 and 5 addresses forecasting CGPP efficiency at different loads and solar irradiance under varying weather conditions. In Section 6, the modeling of various components of the power network and problem formulation with constraints are given. The development of various simulation models suitable for the case study industry with the firefly algorithm and E-Tap software is given in Section 7. Section 8 highlights the simulations of different cases and the results with their brief interpretations. The results and discussions are mentioned in Section 9. Finally, conclusions with certain RES integration guidelines are drawn in Section 10.

Description of Industry under Case Study
In this paper, the differentiation of the proposed heuristic method compared with similar works is established. In support of this work, the case study of a known large-scale continuous process chemical industry; KCR Oxygenated Chemical Complex situated in Andhra Pradesh of India with location latitude of 17.90 N and longitude of 80.80 E. The industry has basically two units of the main process plant which has a big water treatment plant, electrolysis units and the cogeneration power plant. It has both CGPP and SPP and is used for modeling and optimization. This section describes the power distribution network, the CGPP setup, and the details of the solar photovoltaic power plant (SPP).
The process plant requires high-pressure steam of 32 kg/cm 2 at higher temperatures of 240 • C for which a coal-fired steam boiler is installed. The surplus or waste steam is used to run a turbine-generator set to generate power. The total power requirement of the process plant is 36 MW, of which a process load demand of 30 MW is required during normal operations of the industry. The remaining 6 MW is always necessary to support base loads such as water systems, vacuum system and ventilation in plants. When the CGPP is off, the base loads of 6 MW are supplied by the utility grid with the help of a demand management controller, which is used to curb the unwanted demand spike due to the sudden turn-on and turn-off of the large induction motors. Since the industry is continuously operated, the net load of the plant is almost constant. Moreover, some of the loads that run-in daytime are turned off during the nighttime. However, the lighting load offsets the net load during the night. Hence, the load profile of a full day is constant. For reactive power compensation, a shunt capacitor with a rating of 10.6 Mvar is connected to a 6.6 kV bus.
CGPP with operational problems: a CGPP with two steam boilers and turbo generators (TGs) is installed in the considered industry. The steam generator consists of a radiant boiler with integral superheaters, a forced flow section and a tubular air heater with two forced draft fans and two induced draft fans. The turbines are extraction condensing type with turbo generators (TGs). The TGs are 2-pole, 3-phase, air-cooled and connected to the steam turbine directly. Their maximum output is rated at 30 MW (35.294 MVA at a power factor of 0.85) at 6.6 kV. Additionally, the CGPP is connected to a utility grid for startup purposes. It is worth noting that from an operational point of view, the minimum power output of each turbine would be at least 60% of the 30 MW (i.e., 18 MW) [27]. Hence, both TGs together supply 36 MW, which is the exact active power requirement of the industry. Later, the industry installed a 12 MW peak capacity SPP to meet the statutory requirements and to encourage renewable energy (RE). Furthermore, during the SPP generation time (approximately five to six hours a day), the operation of the CGPP is kept running at a lower output of approximately 26 MW (i.e., 13 MW from each TG), and the SPP supplies the rest to satisfy the power balance. Thus, this lowering of the CGPP's output decreases the overall efficiency of the steam turbine from 34% to 32%, and even the extraction of steam into the process system becomes complex. The generation cost of energy is INR 5.6 per kWh. SPP with operational issues: a few years ago, a 12 MW SPP was installed in the industry. The total SPP was divided into two divisions with each 6 MW. A total of 47,060 PV modules with a peak power rating of 255 W; open circuit and nominal voltages of 37.6 V and 30.5 V, respectively; short circuit and nominal currents of 8.95 and 8.42, respectively; a panel efficiency of 15.5% and an annual power derating of 4% [28] were installed with the required inverters. Power evacuation was performed at the 6.6 kV level. The levelized energy cost was estimated to be USD 0.09 or INR 6.5 (INR is Indian rupees in April 2021).
The following operational interlocks that are necessary for the joint management of the CGPP and SPP are presented next. If any TG of the CGPP is tripped, the industry continues to run with support from the utility grid and SPP when it is available. Suppose the process load is reduced suddenly due to any abnormality in process operations. In that case, there arises a situation of islanding of the CGPP during which the SPP is kept 'off' to avoid unauthorized power from being exported to the grid and to prevent a sudden voltage rise at plant buses. For this purpose, an additional mechanical interlock is provided in circuit breakers (CBs). The moment that the CGPP does not generate power, the SPP trips within a maximum time of 0.08 s (i.e., the trip time of CBs at the CGPP and SPP together). This is well within the limits of 0.16 s prescribed by [6]. Even though this type of local PV integration reduces the power losses by 3% to 9% [29], this scenario indicates a certain curtailment of renewable energy, i.e., PV in the present case, for the efficient and secure operation of power systems that are integrated with grids [30]. The design efficiency of the CGPP is 34% up to 16 MW, which decreased to 32% if the load on the CGPP is 14 MW. This is primarily due to the operation of the SPP. With the abovementioned operational problems of the CGPP and SPP, there is a need to determine an optimal solution for such types of industries. It is also necessary to study the forecasting techniques of the irradiance and efficiency of CGPP.

Forecasting the Efficiency of CGPPs
In CGPPs, estimating the overall efficiency is a complex problem for researchers and system engineers since there are many dependent and independent variables such as the steam extraction rate, heat rate and actual system load. The overall efficiency of the CGPP can be obtained either with the heat rate or percentage load as in (1) [31].
The overall efficiency of CGPP = Heat equivalent of electrical output/input energy from fuel combustion where a is a multiplication factor of load and b is a constant.
For CGPPs, no standard function of the relationship between the percentage load and overall efficiency is available except either performance curves of the steam generator and turbines in relation to the heat rate, temperature, etc., by manufacturers or approximate estimations by researchers. New techniques are developed by the researchers for lower coal consumption in steam boilers but it was observed that the coal consumption rate is too low at lower loads and suggested for further studies [32]. Hence, there is a need to develop a function for this, which can help field operators conduct better load management. The field data show that there is the possibility of obtaining a relation function based on an independent variable, with the help of regression analysis tool box in matlab. The statistical and Machine Learning Toolstrip facilitates nonparametric regression for prediction of responses to the new variables with the help of a trained model.

Terms Used in Regression
Some of the metrics used in regression analysis are briefly explained here, with the equations taken from the standard literature.

a.
MAE: mean absolute error (MAE) is the absolute difference between the true values and the predicted values, as in (1). Since it is absolute, any negative sign in the result is to be ignored. It takes the average of error from each sample.

MAE = True values − Predicted values
where, y is true value,ŷ i is predicted value and n is number of observations. It is generally used for continuous variable data and the lower value of MAE gives a better regression model. b.
MSE: mean squared error (MSE) is estimated by considering the averages of the squares of the difference between the true values and the predicted values of the dataset, as in (2).
Here by squaring the difference between the true values and the predicted values, the higher error value can be penalized. It is the most widely used metric when the dataset contains too high or too low values, but not useful when the data contain huge noise. c.
RMSE: root mean squared error (RMSE) is the square root of MSE, which gives better accuracy of the regression model, as given in (3).
It is more useful when large errors are present in the dataset. It is also the standard deviation of errors. d.
R-Squared: it is the co-efficient of determination and indicates the qualitative information of a dataset and is shown in (4).
where, SSEM is the sum of squared errors by mean line and SSER is the sum of the squared errors by regression line. Here the amount of error is reduced since a regression model is used instead of the mean line. Its value lies between 0 and 1 where the value 1 or close to 1 gives a better performance regression model, i.e., a perfect model and a value equal to zero indicates that the model does not fit for the given data. e.
p-value: it indicates the statistical significance of coefficients relationship of a model. A lower valuer value of less than 0.05, i.e., 5% indicates a meaningful addition to the dataset. It tests the null hypothesis of each term of coefficients. A lower p-value means a significant change in prediction value for a change in response variable. It gives a direction of which term to be retained in a regression model. f.
Intercept: it is a constant and it is the value of dependent variable when all the independent values are zero. g.
Standard error (SE): it is the estimate of standard deviation (SD), as in (5).
where σ is the SD and n is the number of samples. h. t stat: t statistic is the coefficient divided by its SE and is preferred when the sample size is small. It is used in statistics to take a decision of rejecting a null hypothesis.

CGPP Efficiency Estimation
For estimation of CGPP efficiency with respect to percentage load on it, the real-time field data method is used. The available data points of the % load and % overall efficiency is shown in Table 1, which also includes the prediction results. These details are used for training regression analysis to predict the overall efficiency of the CGPP and to obtain the relation function. The linear regression-linear model results are as follows: MAE is 0.2018, i.e., lowest, MSE is 0.2313, R-squared is 0.98 which is higher, MSE is 0.05351, and p-value is 5.2336 × 10 −6 for the load and 4.1566 × 10 −9 for the intercept which are very low. All the regression results are shown in Figure 1. The prediction plot, along with derived plot that shows the intercept at zero loads, is shown in Figure 2.
The coefficients obtained in the linear regression analysis are used for predictions, and the related function for the overall efficiency of the CGPP is shown in (6).
The linear regression-linear model results of the RMSE of 0.2313, R-squared of 0.98, MSE of 0.05351, MAE of 0.2018 and p-values of 5.2336 × 10 −6 for the load and 4.1566 × 10 −9 for the intercept indicate that the regression model is perfectly fit for this problem. The obtained function is validated with the available real-time field data of the CGPP and confirms that the relative error is less than 1%, as shown in Table 1, which is a highly encouraging result. Hence, this function is more accurate and can be used for the CGPP in further analysis and to design suitable controllers for its optimum operations, up to the lowest acceptable load of 50%. For operations below 50%, most CGPPs cannot operate efficiently and are not advisable.
where σ is the SD and n is the number of samples.
h. t stat: t statistic is the coefficient divided by its SE and is preferred when the sample size is small. It is used in statistics to take a decision of rejecting a null hypothesis.

CGPP Efficiency Estimation
For estimation of CGPP efficiency with respect to percentage load on it, the real-time field data method is used. The available data points of the % load and % overall efficiency is shown in Table 1, which also includes the prediction results. These details are used for training regression analysis to predict the overall efficiency of the CGPP and to obtain the relation function. The linear regression-linear model results are as follows: MAE is 0.2018 i.e., lowest, MSE is 0.2313, R-squared is 0.98 which is higher, MSE is 0.05351, and p-value is 5.2336 × 10 −6 for the load and 4.1566 × 10 −9 for the intercept which are very low. All the regression results are shown in Figure 1. The prediction plot, along with derived plot that shows the intercept at zero loads, is shown in Figure 2.    The coefficients obtained in the linear regression analysis are used for predi and the related function for the overall efficiency of the CGPP is shown in (6).

Forecasting of Solar Irradiance with Different Weather Conditions
For satisfactory operation of electrical power systems with DG, optimization is generally performed using prior information on the load profile and weather forecasting data. The prediction of SPV generation either an hour or a day ahead from weather forecasts can greatly help power system operators plan their generation and dispatch. The prediction should be automatic and more accurate by considering seasonally varying weather conditions, which is possible with the help of the machine learning technique used by many researchers. Although, beta distribution is used for solar irradiance modeling, authors in [33] used Weibull distribution to compute hourly PV power output. However, every model has limitations, and some deviations exist. In the present case study, the PV data is not validated with any of distributions, as pointed out in literature since the data is readily available for machine learning training purposes. The dataset of 122 available points, both from the field and online, for the case study location (latitude of 17.90 N and longitude of 80.80 E) is shown in Figure 3a-f which have mixed correlation and inter dependent with GHI. The ambient temperature, relative humidity and precipitation have a greater impact on GHI, mostly direct proportional up to 40 • C and starts falling beyond this point. The dew point is directly proportional to GHI, whereas the wind speed has some inverse relationship and mostly no concluding relation to the sky cover. Hence it can be concluded that the correlation between GHI and other parameters has a mixed relation. Some of the sample values are shown in Table 2. These details are used for training regression analysis to predict the solar irradiance, i.e., the dependent variable.
The results of the trained model are shown in Table 3.
The coefficients obtained show that the 'p' value of sky cover is 0.601, which is higher than the acceptable limit of 0.05 [34], i.e., it is not significant at the 5% significance level and indicates that this parameter does not have a significant effect on irradiance. Hence, the regression model is retrained without sky cover. The results are shown in Table 4, and the prediction plot is shown in Figure 4. In this final regression model, the 'p' values of the parameters are less than 0.05 and indicate that the selected model is suitable for the location under study and may require adding sky cover to other locations based on the field data and geographical conditions. pendent with GHI. The ambient temperature, relative humidity and precipitation have a greater impact on GHI, mostly direct proportional up to 40 °C and starts falling beyond this point. The dew point is directly proportional to GHI, whereas the wind speed has some inverse relationship and mostly no concluding relation to the sky cover. Hence it can be concluded that the correlation between GHI and other parameters has a mixed relation. Some of the sample values are shown in Table 2. These details are used for training regression analysis to predict the solar irradiance, i.e., the dependent variable.   The results of the trained model are shown in Table 3.    The coefficients obtained show that the 'p' value of sky cover is 0.6 than the acceptable limit of 0.05 [34], i.e., it is not significant at the 5% and indicates that this parameter does not have a significant effect on the regression model is retrained without sky cover. The results are sho the prediction plot is shown in Figure 4. In this final regression mod the parameters are less than 0.05 and indicate that the selected mode location under study and may require adding sky cover to other loca field data and geographical conditions.
The regression results of the RMSE of 38.797, R-squared of 0.97, MSE of 1505.2, MAE of 28.749 and p-value of less than 0.05 indicate that the regression model is perfect. The obtained function is validated with real-time field details of the SPP and confirms that the relative error is between 16% and −26%, with the majority of values being less than 5%. The results also confirm a mixed relationship between GHI and other varying parameters. However, there is no confirmative relationship between GHI and sky cover, which is also confirmed from the Figure 3. Hence, this function is more accurate and can be used for solar irradiance prediction in the field, especially for integrated power systems with CGPPs. This can help CGPPs be prepared to maintain the stability of turbo generators under varying weather conditions. Figure 5 depicts the power schematic diagram or One Line View (OLV) of the industry whose operational philosophy is described in subsequent sections. The following guidelines are considered in modeling industrial power networks with PV generation. Figure 5 depicts the power schematic diagram or One Line View (OLV try whose operational philosophy is described in subsequent sections. guidelines are considered in modeling industrial power networks with PV

•
In the OLV and E-Tap, even though a single-phase system is shown networks, the calculations and simulations are performed for a three only.

•
The branch flow equations of the distribution network are formulate of LinDistFlow equations proposed in [35] since the case study distrib is radial.  The modeling of the various components of the system is shown in ( equations are developed based on the following conditions and objectives In CGPP, the active output power of each TG shall not be less than 16 efficiency and not be more than the plant maximum load of 36 MW, un power export agreement with either utility grid or other neighboring indus dition shall also be followed by respective reactive power control for a bett and voltage profile with the help of available capacitor banks. Additionally has certain limits of its maximum capacity beyond which it cannot genera for any additions to the existing PV plant, a further power export agreem as stated above. i. CGPP:

•
In the OLV and E-Tap, even though a single-phase system is shown for three-phase networks, the calculations and simulations are performed for a three-phase system only.

•
The branch flow equations of the distribution network are formulated with the help of LinDistFlow equations proposed in [35]  The modeling of the various components of the system is shown in (4) to (7). These equations are developed based on the following conditions and objectives of the industry: In CGPP, the active output power of each TG shall not be less than 16 MW for better efficiency and not be more than the plant maximum load of 36 MW, unless there is a power export agreement with either utility grid or other neighboring industries. This condition shall also be followed by respective reactive power control for a better power factor and voltage profile with the help of available capacitor banks. Additionally, the PV power has certain limits of its maximum capacity beyond which it cannot generate. In addition, for any additions to the existing PV plant, a further power export agreement is required as stated above. i. CGPP: Here, k CGPP is 5.60 INR.
Here, Equation (8) describes the CGPP active power limits, Equation (9) formulates the reactive power limits of CGPPs for time t, and (10) shows the total cost of CGPP energy.
ii. SPP: Here, k PV is 6.50 INR per kWh.
Here, Equation (11) describes the SPP active power limits, Equation (12) formulates the active power limits of SPP for time t, and (13) shows the total cost of SPP energy.
iii. Distribution Network Modeling: Here, k UG is 7.20 INR. Here, Equations (14) and (15) models the net UG power availability in which power is either exported when there is a surplus during SPP peak generation or imported based on plant demand. Equation (16) captures the net voltage at bus 2 (process load bus) after allowing the line impedance drop. The voltage at bus 1 (UG input) is fixed at 1.0 p.u. since there are no users here (17). Equations (18) and (19) describes the UG active and reactive power limits, respectively. Additionally, Equation (23) captures the net voltage at bus 3 (base load bus) after allowing the respective line drop. Equations (24) and (25) stipulates the voltage limits.
Boundary conditions: Since there is no power flow going out of bus 3 (base load bus), the active and reactive power going out from this bus are zero, as in (26).
In the OLV, P2 and P3 represent the net power at the respective node after considering the respective branch power loss. Various technical parameters and their bounds considered in the case study industry are mentioned in Table 5. As shown in the table, the CGPP runs with satisfactory efficiency with loads from 32 to 36 MW. The constraints of this system extend 30 MW to process loads, i.e., process plant running conditions; an essential load, i.e., a base load of 6 MW irrespective of the plant running status. Additionally, a voltage of 100% is considered at bus 1, i.e., v1, since there are no loads on this bus. However, the voltage limits specified in the IEEE Standard are applicable for v2 and v3 at buses 2 and 3, respectively, since loads are directly connected to these buses. The optimal power flow with the lowest cost solution fulfilling the various constraints above, such as the minimum and maximum power flow limits from the CGPP and SPP, voltage limits and total operating costs, for the case study industry have to be tested. The optimal power flow problem (OPF) (P1) of a cogeneration plant integrated with PV for the minimum total energy cost is formulated as in (27), subject to satisfying the parameters in (8)- (26).
Here, the main objective is to reduce the costs of CGPPs, PVs and UGs subject to the conditions formulated in (8)- (26). Hence, the goal is to check the availability of power over time and to utilize the lowest priced energy at that moment. Here, Problem P1 is linear since the objective and constraints are linear. Although the efficiency of the CGPP is not considered in optimization, since it is linear, it can be easily included in the optimization problem. We conducted a post-efficiency analysis after solving the problem. Moreover, since the case study industry has no storage battery system, it is excluded from the optimization. The details of real-time operating conditions with constraints, decision variables and objective functions of the case study industry are mentioned in subsequent sections.

Optimal Power Flow of PV Integrated CGPP Power System
To run the power distribution network efficiently with cost-effectiveness, a suitable optimization technique is required. There are many optimization techniques such as ant colony and artificial bee colony optimization for electrical distribution works. These methods are developed with inspiration from biological evolution and use the mechanisms of (i) reproduction, (ii) mutation, (iii) recombination and (iv) selection. The ant colony technique is inspired by the foraging behavior of ants and is being used in large optimization problems. Artificial bee colony optimization is being used for continuous optimization work. Particle swarm optimization (PSO), which is a population-based optimization, suffers from the disadvantages of easily falling into local optima with a lower convergence rate. Concerning the widely used genetic algorithm (GA), even though there are some advantages, such as complex problems and parallelism, they suffer from disadvantages, such as the requirement of a careful selection of a new population and more time consumption. In this work, optimization techniques such as the firefly algorithm (FA), the modified firefly algorithm (MFA) and the Excel solver are used. These methods are described in subsequent sections. Additionally, for simulations in this paper, a laptop with an x64-based 11th generation Intel core i5-1135G7 at 2.40 GHz CPU with 16.0 GB RAM and a 64-bit operating system is used with MS Office 365-2019 and MATLAB 2018a. i.
Firefly algorithm (FA): this algorithm was developed by Xin-She Yang in 2008 [36] and is inspired by the flashing patterns and movement of fireflies at night. There are three basic rules in the FA: (i) fireflies are attracted to others irrespective of their sex, i.e., unisex; (ii) the firefly's brightness is determined by the value of the objective function; (iii) the attractiveness is proportional to the brightness, which is inversely proportional to the distance between them. Therefore, a dimmer firefly always moves toward a brighter firefly and moves randomly if there is no bright firefly. The FA is becoming more popularly used due to its population subgrouping capability for local attraction. The FA has been tested in the proposed code for the design of a compression spring to determine the spring's suitability in engineering applications [37]. The FA is governed by Equations (28) and (29) [38].
Light intensity, I = I 0 e −γr 2 (28) Firefly attractiveness, β = β 0 e −γr 2 (29) Here, I 0 is the initial brightness value at the source, i.e., distance r = 0; β is the original light attractiveness between two fireflies at i and j, i.e., at r ij = 0. The Cartesian distance is used to estimate the distance between any two fireflies i and j (30).
where 'd' is the number of dimensions, and x i,k is the kth component of firefly i. The updated location of firefly i to a brighter firefly j is given by (31) [39], where α is the scaling parameter, which controls the step size, and rand is a random number drawn from different distributions, such as Levy flights and the uniform distribution. In the above equation, the first term indicates the present location of the firefly, the second term controls the attractiveness and the third term describes the random movement if there are no bright fireflies. We used penalized methods to address the constraints as in (32) and (33) and denoted as x in (34). where, The matrices A eq and B eq , u b and l b are described in Appendix A. Bound projection scheme: If the computed x from the algorithm is lower than l b , then we project x to be equal to the lower band; conversely, if the computed X from the algorithm is greater than u b , then we project x to be equal to the upper band. ii. Modified firefly algorithm (MFA): even though the FA is preferred in most microgrid applications due to its simplicity and good efficiency, it suffers from certain drawbacks such as poor exploration and becoming trapped in the local minima. This calls for an improvement with hybridization, which has been addressed by many researchers [40,41]. In this paper, the FA is modified by augmenting it with the Levy flight distribution and penalty approaches. Levy flight (LF): LF was introduced by Paul Levy in 1937 for statistics [42] and is now being used in economic, physical, biological and engineering fields. Levy motion with characteristics of a long tail or many small steps (i.e., flights) and frequent jumps leads to new clusters, i.e., generating new generations at distances that are randomly distributed. The new generations are closer to the previous optimum population, i.e., speeding up exploitation. Finally, the new generation is evaluated for its best fit. An example of 1000 steps of an LF with (0,0) as the starting point is shown in Figure 6. The direction is uniformly distributed with a Levy distributed step size having values of α = 1 and β = 0. Levy flight can be used to control firefly movement. The random step length, drawn from a Levy distribution, is shown in (35) [43].
where u and v are normally distributed as in (36) and (37).
iii. HyPROS: This add-in MS Excel solver, having a what-if analysis facility, is being used to determine either the maximum or minimum of a formula in a cell, subject to its constraints. In this solver, three solution methods are available: (a) nonlinear GRG (generalized reduced gradient), which is very fast and mostly used for smooth equations to obtain a local solution in milliseconds; (b) LP (linear programming) to solve global problems in moderate time; (c) the evolutionary method, which is based on natural selection and is used to acquire a global solution of non-smooth and nonlinear functions, which takes more time. In LP methods such as LinProg of MATLAB or Excel solver's GRG and LP, the main objective is to minimize a linear objective function over certain decision variables and linear constraints. However, there is no guarantee of obtaining a global solution with LP techniques, but the Excel evolutionary method can give a good or adequate solution. Hence, researchers seek other options for improved methods.
If X < lb, X = lb. Otherwise, if X > ub, X = ub end If the computed x from the algorithm is lower than lb, then we equal to the lower band; conversely, if the computed X from the algor than ub, then we project x to be equal to the upper band. ii.
Modified firefly algorithm (MFA): even though the FA is preferred in m applications due to its simplicity and good efficiency, it suffers from backs such as poor exploration and becoming trapped in the local min for an improvement with hybridization, which has been addressed searchers [40,41]. In this paper, the FA is modified by augmenting it flight distribution and penalty approaches. Levy flight (LF): LF was introduced by Paul Levy in 1937 for sta is now being used in economic, physical, biological and engineering fi tion with characteristics of a long tail or many small steps (i.e., flights jumps leads to new clusters, i.e., generating new generations at dis randomly distributed. The new generations are closer to the previous ulation, i.e., speeding up exploitation. Finally, the new generation is ev best fit. An example of 1000 steps of an LF with (0,0) as the starting po Figure 6.

Numerical Simulations and Results
The case study of industry with its OLV and technical parameters given in the preceding sections is simulated with all the above referred methods for the objective function for two irradiance values of nil and a maximum level of 100% (in the FA and MFA) and for one in the solver.
Case 1-firefly method: here, the classical FA is implemented in the MATLAB software for the case study industry for a 24 h day with a time gradient of one hour. The simulation results are shown in Table 6. These results are in line with the objective function of the industry concerning the power balance and meeting the voltage limits prescribed by IEEE 1547-2020. Additionally, the obtained results are better than existing industry practices. Table 6. FA results (-ve sign indicates the export of power). Case 2-modified firefly method: here, the MFA is implemented in the MATLAB software for the case study industry for a 24 h day with a time gradient of one hour, the same as that of the FA simulations. The simulation results are shown in Table 7, and the fitness curve is shown in Figure 7, which indicates perfect convergence reaching the optimum solution within 100 iterations. Case 3-solver: for the simulation, all three models of the solver method are used with the maximum PV contribution. A summary of the results obtained with the MS Excel solver, i.e., HypROS, for a particular time period of one hour, where the solar irradiance was maximum, is shown in Table 8.   Table 7. MFA results (-ve sign indicates the export of power). Case 3-solver: for the simulation, all three models of the solver method are used with the maximum PV contribution. A summary of the results obtained with the MS Excel solver, i.e., HypROS, for a particular time period of one hour, where the solar irradiance was maximum, is shown in Table 8. To get the values of the real power in watts, the p.u. results of Table 8 were multiplied with base MVA which is 40 here. Therefore, the results with nonlinear GRG and Simplex methods for active power of UG is 8.4 MW (export to UG), CGPP is 36 MW and PV is 8.4 MW. Whereas the active power of UG is 6.4 MW (export to UG), 34 MW from CGPP and 7.2 MW from PV plants with the evolutionary method; here the losses are more due to reduced voltage. The voltages given in the table are multiplied with base kV of 6.6 to get the real values; which is 6.798 kV with nonlinear GRG and Simplex methods and is 6.534 kV with evolutionary method. Additionally, the time taken for simulation is 0.02, 0.03 and 1.2 s for nonlinear GRG, Simplex methods and nonlinear GRG and evolutionary methods, respectively.

Hour P UG (MW) Q UG (Mvar) P CGPP (MW) P PV (MW) Q CGPP (Mvar) v 2 (kV) v 3 (kV)
The results of the GRG and LP methods are the same whereas the evolutionary method gives different results. The results of GRG and LP are very encouraging, and this technique, mostly with LP, can be used by industries that do not have expertise in optimization techniques. Hence, the method has been converted into a ready-reckoner type solver under the name HyPROS (Hybrid Power Resources Optimization Solver) and shown in Appendix B. The GRG or LP methods, which have given faster responses, can be adopted by industries. In this paper, the words solver and HyPROS are interchangeably used.

Results and Discussions
The obtained results are discussed and tested for some of the known key performance indicators as follows: i.
Power sharing and voltage profile: the results obtained for all the techniques used in this work are summarized in Table 9. Here, sample PV penetration of zero and a peak level of 100% are considered for comparison purposes. Additionally, the active power sharing between different resources by the FA and MFA, along with the bus 2 voltage, is shown in Figure 8. The active power sharing curves show that the power allocation and voltage profile is smooth without any rapid swells and sags in MFA due to Levy flight when compared with the FA. Hence, the local CGPP cannot be affected much and exhibits stable operation. The MFA curve is better than the FA curve. Additionally, the voltage profile is comparatively better in the MFA model. It shows that all the results are almost the same with a better voltage profile with the MFA model. Additionally, the Excel solver results are equally comparable with other optimization techniques and can be used in the field. The energy cost per hour of the existing practice is INR 0.2304 million. However, the cost estimates are reduced by 9-15% with various optimization techniques. Even though the best estimates are obtained with the solver technique, the next-best MFA is the most useful for the field due to its technical advantage of a better voltage profile. This saves approximately INR 228 million, i.e., 11% savings annually with the MFA optimization model of the energy mix. ii. PV penetration ratio: here, we estimated the PV penetration ratio for a 6.6 kV system to check for the penetration levels, as in (38).
PV Penetration ratio = S PV f eeder/n loads S peak (38) where S PV feeder is the PV power installed for n loads with the estimated peak power of S peak . In the present case study, the PV penetration ratio = 100 (12 MW)/(1*36) = 33.33%, which is termed 'high penetration' [44,45]. However, from past field data, the maximum PV power output was 8.5 MW, i.e., at an average performance ratio of 70%. In such cases, the PV penetration ratio is adjusted to 23.6%, which is marginal. iii. Amount of energy imbalance (AEI): the authors in [46] proposed the AEI as in (39), which represents financial amounts but can be used for power systems and indicates the difference between the power supplied by and imported from grids.
where SA is the settlement amount (power produced), SB is the settlement basis withdrawal amount, AP denotes the amount of purchase with a bilateral agreement (BA), BDAM is the buying from the day-ahead market (DAM), SDAM is the amount of sales to the day-ahead market, ALS is the amount of load shedding, and AUL is the amount uploaded. In the present case study, for the time period of 10 to 11 AM  (Table 7), the PV generation was 6500 kWh with a load of 36,000 kWh, CGPP generation was 32,000 kWh and export to UG was 2510 kWh.
AEI = (36,000 + 6500) − 0) + (0 − 0) + (0 − 2510) + (0 − 36,000) = 3990 Here, the positive AEI indicated a power surplus and, hence, exporting. In the case of a negative AEI, there will be a shortage of power and imports from the grid. This calculation can help develop better control strategies since we cannot destroy the energy. iv. Validation with E-Tap: LinDistflow linearizes the power flow equations [47], and the computed 'X' may not be feasible for actual nonlinear solvers, e.g., E-TAP. Therefore, the results obtained in HyPROS, from FA and MFA are validated with E-Tap, the power system analysis software [48], to determine the voltage performance at the optimum power flow conditions. Here the optimum power flow conditions obtained in those algorithms are used to simulate the power flow conditions in E-TAP. On simulating these optimum power flow values in E-TAP, voltage and power factor at various buses are checked for their validity in the system. The response from E-Tap for v 2 and v 3 confirms that of the various optimization method results are well within the limits prescribed in IEEE 1547-2018. A summary of the validation results is shown in Table 10. An illustration of the E-Tap load flow result for the nonlinear GRG solver is shown in Figure 9. Here, the single line diagram or OLV of the power network is transformed into an E-Tap diagram. E-Tap estimations confirmed that the generator was operating with the best power factor of 92.74%. v.
Voltage Ramp Index: it is a known fact that with the increased PV penetration, the voltage stability of grid will be affected [49]. In addition to the abovementioned performance parameters, we propose the 'voltage ramp index (VRI)' as in (40), to determine the hourly or any step-time voltage variations or spikes in a power distribution system due to varying PV power output. In the selection of the maximum value, the sign of either +ve or -ve should be ignored to obtain the index value. (40) where v t+1 , v t , v n+1 and v n are the selected node voltages at t + 1 th , t th , n + 1 th and n th hours, respectively. For explanation purposes, the FA results in Table 6 are reproduced here in Table 11 with hourly voltage readings at bus 2. The step voltage or ramp voltage is the change in voltage from the preceding hour to the present hour, e.g., for the 14th hour, the step or ramp is v 2 14 − v 2 13 = 6.63 − 6.71 = −0.08 kV per hour, which indicates a voltage drop of 0.08 kV. Therefore, for the selected 24 h day, the VRI is 0.12, which is due to the PV power being 'off' at 6 PM, i.e., at sunset. Such information is useful in selecting voltage regulators or load tap changer adjustments. Additionally, the MFA results in Table 7 are reproduced here in Table 12 with hourly voltage readings at bus 2. The step voltage or ramp voltage is the change in voltage from the preceding hour to the present hour, e.g., for the 14th hour, the step or ramp is v 2 14 − v 2 13 = 6.76 − 6.76 = 0.00 kV per hour, which indicates no change in voltage.
For the selected MFA Therefore, for the selected 24 h day, the VRI is 0.7, which is due to the PV power being 'off' at 6 PM, i.e., at sunset. With MFA, it is seen that the number of zero voltage steps are 13 against two in FA, which indicates a smooth voltage profile in MFA and this model can be adopted for industries having cogeneration power plant integrated with large size PV power plants.       v. Voltage Ramp Index: it is a known fact that with the increased PV penetration, t age stability of grid will be affected [49]. In addition to the abovementioned mance parameters, we propose the 'voltage ramp index (VRI)' as in (40), to de the hourly or any step-time voltage variations or spikes in a power distribut tem due to varying PV power output. In the selection of the maximum value, of either +ve or -ve should be ignored to obtain the index value. where vt+1, vt, vn+1 and vn are the selected node voltages at t + 1 th , t th , n + 1 th hours, respectively. For explanation purposes, the FA results in Table 6 ar duced here in Table 11 with hourly voltage readings at bus 2. The step vo ramp voltage is the change in voltage from the preceding hour to the prese e.g., for the 14th hour, the step or ramp is − = 6.63 − 6.71 = −0.08 kV p which indicates a voltage drop of 0.08 kV.

Conclusions
The optimal operation of a grid-connected cogeneration power plant integrated with a solar photovoltaic power plant was determined. For this work, a case study of a known continuous process chemical industry of a similar type was selected. Forecasting techniques for cogeneration efficiency versus its operating load were proposed and obtained a satisfactory relation between its efficiency and load, which is useful for cogeneration operators to run it more efficiently. Additionally, solar irradiance forecasting techniques were proposed and obtained a reasonable relationship between the irradiance and varying weather parameters such as ambient temperature, wind speed, relative humidity, dew point, etc. They were verified real time field results and found to be encouraging. This technique will be useful for field operators to forecast the solar irradiance and in turn solar power, for any type of seasonal weather conditions combination. To solve the optimization problem of various energy resources, different optimization techniques, such as the solver, firefly algorithm and modified firefly algorithm, were proposed and implemented. The simulation results of the MFA method were compared with those of the FA method and confirm that the MFA method is more effective and well suited for industries with grid-connected cogeneration power plants along with locally integrated solar PV power systems. The set points which were computed from the optimization were applied at the field and obtained a satisfactory CGPP efficiency of above 34%, which is the lowest limitation given by the steam boiler manufacturer. Additionally, a handy optimization solver was developed for industrial use. This can be used by any technical person, not having knowledge of complex optimization techniques, who can easily use the solver which basically requires experience in MS Excel. The authors concluded that the MFA method, in addition to the solver, is more useful for similar types of industries and maximizes the harnessing of solar energy.
The results obtained with various optimization techniques were validated with E-Tap, a standard power system analysis software, for voltage profile and were found to be satisfying the voltage limits prescribed in the latest IEEE standard 1547. With this satisfactory technical validation, the same can be comfortably implemented at the field level. Various performance indicators were also used to confirm the suitability of the techniques adopted and obtained encouraging confirmations. Additionally, a new performance indicator, voltage ramp index was proposed to know the voltage steps due to a varying mix of energy resources. The voltage profile obtained with FA and MFA optimization techniques was verified with the newly proposed voltage ramp index and useful information was obtained. This work mainly focused on the optimal operation of PV integrated cogeneration plants and integration problems. The other major integration problems of reactive power and protection systems were addressed by the authors separately. Initially, there was some comprehension that the higher penetration of RERs can disturb existing electrical networks. However, continuous research on RE integration problems and advancement of artificial intelligence techniques into the electrical power sector can comfortably address these problems. In future works, nonlinear forecasting methods for a year of PV systems and other improved versions of optimization techniques can be considered. The available technologies, such as AI for forecasting demand-supply gaps and optimal utilization of resources and peer-to-peer power trading with blockchain technology, are to be used by industries for the smooth operation of electrical power systems.