Classification and Analysis of Optimization Techniques for Integrated Energy Systems Utilizing Renewable Energy Sources: A Review for CHP and CCHP Systems

Energy generation and its utilization is bound to increase in the following years resulting in accelerating depletion of fossil fuels, and consequently, undeniable damages to our environment. Over the past decade, despite significant efforts in renewable energy realization and developments for electricity generation, carbon dioxide emissions have been increasing rapidly. This is due to the fact that there is a need to go beyond the power sector and target energy generation in an integrated manner. In this regard, energy systems integration is a concept that looks into how different energy systems, or forms, can connect together in order to provide value for consumers and producers. Cogeneration and trigeneration are the two most well established technologies that are capable of producing two or three different forms of energy simultaneously within a single system. Integrated energy systems make for a very strong proposition since it results in energy saving, fuel diversification, and supply of cleaner energy. Optimization of such systems can be carried out using several techniques with regards to different objective functions. In this study, a variety of optimization methods that provides the possibility of performance improvements, with or without presence of constraints, are demonstrated, pinpointing the characteristics of each method along with detailed statistical reports. In this context, optimization techniques are classified into two primary groups including unconstrained optimization and constrained optimization techniques. Further, the potential applications of evolutionary computing in optimization of Integrated Energy Systems (IESs), particularly Combined Heat and Power (CHP) and Combined Cooling, Heating, and Power (CCHP), utilizing renewable energy sources are grasped and reviewed thoroughly. It was illustrated that the employment of classical optimization methods is fading out, replacing with evolutionary computing techniques. Amongst modern heuristic algorithms, each method has contributed more to a certain application; while the Genetic Algorithm (GA) was favored for thermoeconomic optimization, Particle Swarm Optimization (PSO) was mostly applied for economic improvements. Given the mathematical nature and constraint satisfaction property of Mixed-Integer Linear Programming (MILP), this method is gaining prominence for scheduling applications in energy systems.


Introduction
In today's world, the growing realization of climate change has been the driving force behind the increasing trend of renewables share in the energy sector [1]. Renewable energies development has been growing at a fast pace, in which many countries around the world are establishing their future energy scenarios based on such sources of energy [2]. Over the past decade, renewable energy sources including hydropower, wind, solar, and biomass experienced an 8% increase in the final electricity generation, reaching at 28% in 2020 [3]. Solar and wind energies were the primary contributors to this expansion, in which solar PV and wind power experienced 16% and 12% annual growth in the renewable electricity generation share, respectively. The efforts to increase renewable generation resulted in reduction of global carbon intensity, and consequently, CO2 emissions in the power sector by 2.5% and 1.3% in 2019, respectively. In the light of the global pandemic (Covid- 19), energy-related CO2 emissions were estimated to decrease by 8% in 2020, since the global energy demand was projected to drop by approximately 6% [4]. Although the demand for renewables increased by 1.5% in the first quarter of 2020, renewable energy was not the only contributor to the lower emissions. In fact, milder than average weather conditions and lower energy demand resulted in significant lower heating demand, particularly in the Northern Hemisphere region. This reduction was stronger in the industrial sector at 4.2% than in buildings at 1.8% [5]. Realizing this strong impact can illustrate the dependency of energy demand, and consequently, CO2 emissions on the heating demand. Despite such strong dependency, the heating sector has received nominal contributions to reduce carbon intensity. In the heating sector, fossil fuels are the dominant source of energy, being responsible for 73% of the total share. Given that heat dominates half of the energy end-use share, it is of the utmost importance to reduce the fossil fuels intensity in order to step closer towards the net zero-carbon society.
In the context of efficient energy utilization and effective energy management, integration of thermal energy production in power generating systems introduced a promising approach to energy-related issues. This method of energy generation was soon recognized as Combined Heat and Power (CHP) and has been appreciated by many economically developed countries. In CHP systems, the primary objective is improving the efficiency and economic by generating power and heat simultaneously. Achieving this common goal allows such systems to lower the primary energy consumption and environmental consequences. Despite the long existence of cogeneration and trigeneration systems, the introduction of centralized power stations with significantly lower costs turned away the attentions and interests from implementation of such configurations. To this end, developments of CHP and CCHP experienced a significant reduction followed by an increasing trend over the past decade. In the developed countries such as the UK and Canada, imposing financial incentives by the government have created a competitive market for CHP and CCHP systems, resulting in an increasing interest particularly in the industrial sector. Nonetheless, the capital-intensive nature and complexity of such systems still outweighs the advantages offered. To this end, development of CHP and CCHP systems has been mostly considered for the decentralized energy generation in the industrial sector. In particular, the chemical industry is one of the most energy intensive and the biggest utilizer of CHP and CCHP systems. For instance, the chemical industry accounts for the largest portion of energy demand in the UK, in which only 36% of the demand is consumed as electrical energy while the remaining is spent on processing heat and steam [6]. As of 2016, CHP was developed in 53 chemical sites across the UK, which accounted for a significant energy generated in this sector [7]. In the US, a similar trend can be illustrated as the CHP capacity in chemicals and refining sector was recorded at 20,000 MW and 15,000 MW in 2016, respectively [8]. Although the current trend seems appealing, CHP and CCHP systems represent substantial capital investment, suffer from complexity in the control strategy, and require optimized thermodynamic performance in order to be considered for government subsidies. In this regard, it is required to solve the associated problems with both thermodynamic and economic performance of CHP and CCHP systems. Owing to the fact that these two perspectives are tied in any energy system, the thermoeconomic concept was introduced, and it now represents one of the most important aspects of energy system performance.
Optimization is a powerful tool and a promising solution to any operation-related issues. The ever-increasing complexity of engineering systems, growing demand for accuracy, and searching for optimal and robust designs introduce additional difficulties that could only be accounted for by developing optimization models. In theory, optimization is the process of maximizing or minimizing an objective function, by consistently selecting and/or calculating possible outcomes within a certain set of candidates. Energy system optimizations are widely used to provide insights on a variety of different scenarios, when the energy system is subjected to numerous thermodynamic, economic, and environmental constraints. Complexity of an energy model is heavily relied on the mathematical representation and the energy system that it is based on [2]. Even with highly accurate modeling, high quality data, and vast computational resources, the complexity of energy models seeks for a significant amount of modeler judgment to validate model results [9]. CHP and CCHP systems are no exceptions; in fact, highly complicated given the presence of more than one output product in the system. Over the past years, there have been significant contributions on investigating the technical and economic aspects of heat and cooling generation along with power. These contributions are mainly evolved around optimization of thermoeconomic performance with respect to energy balance, heating and cooling demand availability, pricing and financial incentives, environmental policies, and power and heat scheduling. The optimal solutions to such optimization problems are obtained using a wide range of search-based (search method) and mathematical-based optimization techniques (calculus method). The search method follows a random or systematic pattern (sequence) progressively in order to improve the objective function. On the other hand, the calculus method, or famously known as the gradient method, utilizes first and/or second-order derivatives in order to determine the extremum or minimum of an objective function.
There have been a limited number of studies that reviewed the existing optimization techniques for development of integrated energy systems. In one of the first attempts, exergoeconomic optimization of CHP systems and different cost allocation techniques were presented and reviewed by Abusoglu and Kanoglu in Ref. [10]. The authors classified the thermoeconomic optimization methods into two primary groups namely the algebraic method and calculus method. The presented techniques were compared with respect to one of the most famous CHP optimization problems known as Constrained Generalized Additive Model (CGAM) [11]. In a more comprehensive review, the applications of heuristic optimization algorithms for solving the Combined Heat and Power Economic Dispatch (CHPED) problem was investigated by Nazari-Heris et al. [12]. The authors in Ref. [12] discussed the necessary optimization procedures for the solution of the CHPED problem, taking various operating constraints into account. In Ref. [12], more than 18 heuristic optimization techniques were presented, and their performances were analyzed by evaluating the outcomes of the preliminary research studies. Application of CHP systems in microgrids provided the opportunity for the authors (Khan and Wang) in Ref. [13] to review the preliminary optimization research works with regards to the multiagent system (MAS) technology. Khan and Wang investigated the operation of microgrid networks using MAS with different optimization techniques [13]. For the case of CCHP systems, Cho et al. in Ref. [14] studied the energetic and exergetic analysis of CCHP systems, and realized the existing methods in performing such assessments. The study in Ref. [14] was mainly evolved around model development and optimization of CCHP systems. In an alternative approach, modeling, planning, and energy management of CCHP systems for microgrid application was reviewed by Gu et al. [15]. The study in Ref. [15] elaborated on technical, environmental, and economic optimization of CCHP microgrid networks. In a more general approach, Cui et al. in Ref. [16] evaluated the performance of heuristic algorithms in solving multiobjective optimization problems for energy saving application. Similarly, multiobjective optimization of district energy systems utilizing variable renewable generation was reviewed by Wang et al. [17]. The study in Ref. [17] looked into the concept of energy systems integration in particular, and presented a significant literature of energy performance in such systems. The uncertainties associated with renewable generation urge for performance optimization for both short-term and long-term planning. In this regard, Zakaria et al. in Ref. [18] evaluated the stochastic optimizations in renewable energy applications. The objective of their research (Ref. [18]) was to introduce the most recent advancements in stochastic modeling and optimization of energy systems utilizing variable renewable generation. Implementing the optimization algorithms can be done using different open source or private property tools for both research and practical use. In Ref. [19], maturity of open source tools for developments of various optimization techniques was investigated. The authors (Ref. [19]) evaluated more than 30 open-source tools based on 81 different functions that they could offer.
Future energy systems will no longer depend on just the power sector. With the energy transition in progress, vast majority of future energy scenarios are heavily relied on sector coupling. Although the power sector remains one of the key players, it will not account for the largest portion of the demand, even if the fossil fuel-based heating units and transportation means are replaced with electric-based units [2]. There is a need to go beyond the power sector and integrate other forms of energy demand into the energy systems optimization models. Moreover, the energy transition strategy on the long-term is subjected to a combination of factors that are severely uncertain, including socioeconomic impacts, resource availability, global energy efficiency investment, and technology improvement [9]. The global pandemic has caused significant disruption to the energy sector, changing the dynamic and the future of the clean energy transition. The temporary depletion of energy demand and CO2 emissions in 2020 were the results of economic downturn. The global epidemic has the potential to change the priorities both in the private and public sectors, which may strongly affect the dynamic of the energy sector, and delay the sustainable development scenarios for an uncertain period. With this level of uncertainty, there is a need to integrate the current uncertainties in all the future energy systems optimization models. It is required to first model the optimization problems in a more general and local approach, and further solve for the global optimization problems. In this context, it is of the utmost importance to realize the optimal technique for programming and solving the energy systems optimization models.
Understanding the need for the current study was first realized when no comprehensive literature could be found on constrained and unconstrained optimization techniques. There are nominal contributions on realizing the application of each method in optimization of the CHP and CCHP performance. Despite categorizing optimization into the search method and calculus method, optimization problems can be classified into various categorizes based on several factors including, but not limited to constraints, deterministic nature of variables, and the number of objective functions. In this study, optimization methods are classified based on the availability of constraints in the model. Additionally, the possible applications of each technique in CHP and CCHP energy systems are provided and the upsides and downsides are discussed. Moreover, statistical reports of every constrained methodology are presented based on a systematic search. The current work is structured into four sections. Section 1 overviewed the objectives of the current work and reviewed the preliminary studies on energy system optimization models. Section 2 elaborates on unconstrained optimization techniques, and presents the most well known methods in such a type of optimization. Further in Section 3, heuristic and mathematical methods of constrained optimization are discussed and the principle of operation and the advantages of each method are presented. Finally, Section 4 concludes the study.

Unconstrained Optimization Techniques
The field of unconstrained optimization is vast and there exists various techniques for solving such optimization problems. There is no doubt that most structural optimization problems involve constraints to enforce limitations on the system performance. Although, unconstrained models are important to derive the search direction and optimum when optimization models reach to no constraints stage. Such problems can arise from real variables without any constraints, or they can be derived from reformulation of constrained optimization problems. Unconstrained optimization is gaining prominence in linear and nonlinear structural analysis problems [20]. Majority of unconstrained optimization methods employ line search algorithm to search for the optimal solution. In most of the techniques, the process starts by considering an initial point and a suitable step length, which then the termination criterion gets checked at each iteration until it is satisfied. The following methods, including: Nelder-Mead, Fibonacci, golden section, Hooke and Jeeves', Powell's, gradient descent, and coordinate descent, are the most common line search algorithms for solving unconstrained optimization problems [21].

Nelder-Mead Method
The Nelder-Mead approach is one of the most widely used and well-studied gradient-free algorithms [22]. This technique was first developed in 1965 by Nelder and Mead for unconstrained nonlinear optimization problems [23]. It can be employed for multidimensional unconstrained minimization of the objective function. Owing to its simplicity and low computational efforts, it has been widely used for various applications such as power dispatch, distribution state estimation (DES) in renewable energy systems, and optimization of the wake effect in wind farms [24][25][26][27]. This method of optimization utilizes a geometrical shape known as Simplex to search its domain. In its principle, simplex is shaped as 1 vertices in an n-dimensional space, where the worst vertex in iteration gets to be replaced by reflection, expansion, contraction, and reduction [28]. The process in this algorithm starts by initializing an n-dimensional Simplex (guess vertex) with 1 vetices, which then gets to be reshaped and moved in every iteration towards the optimum value ( Figure 1). In Nelder and Mead's original model, four coefficient values were suggested, namely reflection coefficient ( 1 ), expansion coefficient ( 2), contraction coefficient ( 0.5), and shrinkage coefficient ( 0.5). The search carries out in the optimal premises until the stop criterion is met and the vertex corresponding value is returned as the optimal output. In comparison to other unconstrained search-based techniques, the Nelder-Mead method provides a more flexible approach, in which the size and geometry of the simplex can be altered. J.A. Nelder and R. Mead presented their model in Ref. [23] as a more superior and flexible technique compared with that of the Powell method in Ref. [29]. Nelder-Mead optimization is a downhill simplex method for minimizing any function. The application of this technique in todays' systems might not be viable and satisfactory since its capabilities and functions are limited. Although, such a method is considered as the basis of other new techniques and its employment for minimization of the objective function is still considered in a handful of studies [31]. In particular, the Nelder-Mead algorithm was implemented in a thermoeconomic evaluation of the ORC-CHP system (considering the part-load operation) by Lecompte et al. [32]. In their study (Ref. [32]), the objective function was set to minimize the Specific Investment Cost (SIC) through optimizing the size of components (e.g., condenser and evaporator) and thermodynamic parameters (e.g., mass flow rates of air and the working fluid). On the other hand, authors in Ref. [33] employed this method of optimization to minimize the Levelized Cost of Electricity (LCOE). The proposed system by Datas et al. [33] was comprised of solar PV modules, in which the storage of solar energy in the form of heat for power production on demand was evaluated. The multivariable algorithm in this case was evaluated over a matrix of diverse initials so as to avoid any local minimum. Aforementioned, implementation of the Nelder-Mead method in new optimization techniques is feasible at no constraints stage for minimization of a specific objective function along the process. With this regard, Pierobon et al. in Ref. [34] suggested a Generic Algorithm (GA) for multiobjective optimization of a large-scaled ORC-CHP system. In this optimization approach, the Nelder-Mead method was utilized as a suboptimizer to select the tube and shell geometry for a specified velocity in the tubes and on the shell side. The objective of this suboptimizer was to minimize the marginal error of shell speed. Application of the Nelder-Mead method for finding the minimum value of any objective function, makes this technique useful in any optimization model. In Ref. [35], this method was employed by Rybiński and Mikielewicz for minimizing the sum of squares in design optimization of mini-and microchannel heat exchangers. The authors in Ref. [35] developed a statistical model for the determination of thermal characteristics of such heat exchangers. In their design optimization (in Ref. [35]), the Nelder-Mead method was used to find the exponents of Reynolds and Prandtl numbers so as to minimize the sum of squares associated with correction and thermal resistance coefficients.

Fibonacci Method and Golden Section Method
Fibonacci and Golden Section methods are two of the most well-known sequential search techniques, in which the strategy is to utilize the experience gained in the earlier searches in order to obtain the optimal value of the function [36]. In such a type of numerical search method, the values of objective function are found at different combinations of decision variables, and the optimal solution is obtained through a series of criteria. The Fibonacci search technique is a strategy that acquires the smallest interval of uncertainty, where the optimal value is located. The principle of this method is based on the Fibonacci sequence and it is employed to find the minimum or maximum of a unimodal function in a given interval (Figure 2a). With this regard, this method of optimization can also be considered for non-differential or discontinuous objective functions [37]. Despite its fast convergence and no requirements for derivatives, this search method does not obtain the exact optimum value. In contrast, it only provides the smallest interval (final interval of uncertainty) after N iterations, or when it is within the user-specified tolerance. In addition, it should be guaranteed that the optimal solution lies in the initial interval of uncertainty [38]. The golden section method is an approximation of the Fibonacci algorithm that can be employed, when the number of experiments is not specified beforehand (Figure 2b). The difference between the golden section and Fibonacci technique lies in the step length ( ) specification; while in golden section the reduction of intervals is carried out at a constant ratio (asymptotic approximation), the Fibonacci method takes advantage of the inconstant ratio that is based on two consecutive Fibonacci numbers at each iteration. As a direct consequence, in a given number of steps (lower N) Fibonacci is more accurate [36]. Such sequential methods are efficient and accurate for a small number of subintervals, however, due to linear convergence they require more computational time for higher iterations. Additionally, initialization is required to perform this type of optimization [39]. In practice, golden section is the method that is favored since it is simpler to implement and requires no knowledge of the number of function evaluations beforehand [24]. The reliability of these two methods is the primary reason for their utilization in energy system optimization models [20].
Considering the application of the Fibonacci optimization method in CHP and CCHP systems, there are a limited number of studies employing this technique. This can be due to the fact that such an optimization method is more suited for a limited number of iterations and in larger data sets this method tends to output the same final interval as golden section, which is a simpler approach compared to Fibonacci. In this context, Hagos et al. in Ref. [40] investigated the prospect of bioenergy integrated district heating (DH) networks in Norway using the Fibonacci technique. Inspection was carried out for two existing gasification-based plants in different scenarios. One of the objectives was to estimate the power and heat consumption based on the historical data. This was done by implementing a time series approach, in which the Fibonacci method was employed for logistic curve fitting, and the goodness of fit was considered as its optimization criterion. The application of the Fibonacci algorithm in energy systems is mainly exploited in solar PV configurations, in which Maximum Power Point Tracking (MPPT) is required [41,42]. In this regard, such a technique can be advantageous in the future studies that will elaborate on Integrated Energy Systems (IESs) with solar PV integration. Starting with voltage to derive the power output of the system, the Fibonacci method searches for the MPP in the given range.
In contrast to the Fibonacci technique, application of the golden section search method in CHP or CCHP systems is more widespread. Such a technique has been used in many studies to find the optimal value of a unimodal function, particularly in ORC-CHP systems. For instance, Quoilin et al. in Ref. [43] designed a solar thermal system, comprising of parabolic trough collectors coupled with an ORC for small-scale applications. In their study (Ref. [43]), influence of temperature glide on the overall efficiency was investigated. It was presented that the glide temperature had to be set on its optimum value (corresponding to optimal flow rate of the working fluid) in order to obtain the highest efficiency. To achieve this value, the golden section search algorithm was employed. In another ORC-CHP system, which was also mentioned in the Nelder-Mead method, golden section was used to find the maximum net power output under part-load operation [32]. Reviewing the application of golden section search in trigeneration systems reveals that there are several studies contributing to this subject. In both Refs. [44,45], this method of optimization was employed in the Engineering Equation Solver (EES) for a Photovoltaic-Thermal (PV/T)-integrated CCHP system. While in Ref. [45] exergoeconomic analysis was conducted to evaluate and optimize the performance, authors in Ref. [44] considered thermoecological cost assessment of the PV/T-CCHP plant. An alternative approach was taken in another study, where the golden section was implemented as a suboptimizer for performance analysis of a co-firing biomass/NG-CCHP configuration [46]. In the proposed optimization model, minimizing the marginal error associated with energy balance equation for a given gasification temperature was performed with the golden section search method. This technique could also be used for finding the values of independent variables that can result in the specified or optimum values of the dependent variables. In this context, C.Y. Zheng and X.Q. Zhai used golden section to obtain the optimal power generation unit (PGU) capacity under different operating scenarios [47]. A similar approach was taken to find the values of independent variables that could output the specified power break for designing an optimized micro-Stirling engine in tri-generation applications [48].

Hooke and Jeeves' Method and Powell's Method
The Hooke and Jeeves' method, famously known as pattern search, is a sequential technique that was first developed by Hooke and Jeeves [51]. This searching method takes advantage of two different kinds of moves at each iteration including the exploratory move and pattern move. In the first stage, a series of coordinated decent steps are performed in the premises of the current point ( Figure 3a). This is done by adding an increment to one of the components so as to obtain the best possible point in its neighborhood. The second stage considers the achieved data and it conducts a pattern move along the stepping valleys, in which the search carries out between the line joining the first and the last points in a cycle [52,53]. At each step, the stop criterion is checked and if the exploratory move leads to the optimum point, the search is carried out by performing the pattern move; otherwise, the process is reset. In general, the exploratory move intends to achieve useful information regarding the objective function in the neighborhood of the current base point, while the pattern move identifies the best search direction by using the acquired information [54]. This method of optimization is similar to univariate search; therefore, its convergence speed is quite slow, might fail depending on the initial conditions, and it is very sensitive to the parameters' variations. To overcome the problems with this technique, Powell developed a more comprehensive and efficient algorithm considering the same pattern search strategy as Hooke and Jeeves [29]. Powell's method is an extension of pattern search and it is based on the principle of conjugate directions. This widely used and accurate technique takes advantage of previous solution information in order to create a new search direction (Figure 3b). It discards one of the coordinate directions in favor of the pattern direction [49]. The algorithm creates a set of linearly independent directions and it performs unidirectional searches along each direction connecting the two points.
Expectedly, the Hooke-Jeeves' method was not considered for developing optimization models in many studies. Since the approach cannot guarantee global optimality, majority of studies avoid implementing this technique. Hence, its application was mostly contemplated for local minimization or calibration purposes. In Refs. [55,56], this method was employed to calibrated Stirling-based micro-CHP systems in order to minimize the discrepancies (cost function) between experimental and simulation results by Cacabelos et al. and Bengoetxea et al., respectively. The Hooke-Jeeves' algorithm would select appropriate input-parameter values aiming to achieve cost function minimization. In an alternative approach, this method of optimization was coupled with the Particle Swarm Optimization (PSO) technique to develop a hybrid metaheuristic optimization algorithm [57]. The suggested ecoenvironmental optimization model was implemented for a central solar heating plant. In particular, the Hooke-Jeeves' method was employed to avoid the probability of getting trapped in a local optimum. Surprisingly, Powell's technique has not been a favorable method of optimization for coor tri-generation systems as well. Considering its advantages such as accuracy, fast convergence, and reliability, one might hope to find its employment more often. As mentioned in Section 2.1, J.A. Nelder and R. Mead presented the superiority of their model in comparison to the Powell's method. Given the flexibility and robustness of the Nelder-Mead method, it is no surprise if practical implementation of Nelder-Mead is favored over the Powell's method. Although, the superior performance of the Nelder-Mead method is heavily affected as the dimensionality of the problem increases [24]. In this context, there is no definitive comparison between the two methods, concerning energy system optimization models. Depending on the application, the performance superiority may vary from one method to another [58]. With development of new evolutionary methods, the applications of such direct search techniques are fading or modifying. In a two-part research, an ORC-CCHP system was developed and subsequently its thermoeconomic optimization was conducted using Powell's algorithm in EES software [59,60]. A similar approach was taken in Ref. [61] by Shokati et al. to optimize a Kalinabased ORC-CHP system from an exergoeconomic point of view. Bellos and Tzivanidis evaluated and optimized their proposed ORC configuration with various working fluids using Powell's algorithm [62]. Interestingly, the implementation of this method can cover various applications, for instance, it can be employed to develop a design optimization model for energy production units such as Stirling engine [55], or solve the dispatch problem of any CHP or CCHP system at the local level [63].

Gradient Descent Method and Coordinate Descent Method
The gradient descent (ascent) algorithm is one of the oldest first-order derivative optimization methods that is employed for unconstrained local minimization (or maximization) [65]. The strategy of such a technique is based on the calculus method, using gradient points in the direction of the fastest decrease (or increase) for a differentiable function. The process starts by estimating an initial value for minimum (or maximum) design, and further calculates the steepest descent at that point using the gradient of the function (Figure 4a). This method can take advantage of the fixed step size or line search techniques to search along the gradient line at every step. Line search is a technique for identifying a suitable step length, in which a sequence of candidate values is tested until a certain condition is satisfied. The line search continuously looks for the minimum (or maximum) points along the derivative line until the convergence is reached. The line search technique can be time consuming, while the fixed step size method results in poor convergence. The coordinate descent (or ascent) method is also an iterative approach with the difference that, minimization (or maximization) is achieved through the coordinate selection rule, rather than using the gradient that points to the direction of the steepest descent (gradient descent method). By fixating some components at a time, the minimization is achieved using line search along one specific coordinate direction (Figure 4b). With this regard, such a method is applicable for both differentiable and derivative-free unconstrained objective functions [66]. Both gradient and coordinate methods benefit from simplicity and stability, however, both methods rely heavily on the accuracy of initial point and step size. Although there are various methods to accurately identify step size, but they are time consuming and complicated to implement, which consequently results in low convergence speed [67].
Concerning the applicability of gradient descent and coordinate descent methods, not many studies considered their application in CHP or CCHP systems. In addition to downsides mentioned earlier, implementation of such methods is only feasible at the local level. There exist several modern methods that are capable of handling complexity of todays' energy systems, and satisfying the accuracy required at the global level. With

Stop
Yes recent developments in computer technology and machine learning, there has been increased interest in the coordinate descent method for the problems in which computing the gradient is not feasible [68]. In the case of gradient descent, this method has been the foundation of many complex and non-linear optimization techniques such as quasi-Newton's method [20]. Moreover, various extensions of such a technique have been presented over the past decades such as the fast gradient method and accelerated gradient method, with the intension to faster convergence for convex problems [69]. Nevertheless, a handful of studies took advantage of its simplicity for employment in CHP systems. In Ref. [70], the application of the gradient method was considered in optimal operation and planning analysis of an IES by Li et al. Their study (Ref. [70]) investigated operating modes of compressors in terms of energy flow, and calculated the partial derivatives of each state variable using the gradient descent iterative technique. An alternative approach was taken by researchers in Spain, where this method of optimization was implemented to calculate the values of 12 decision variables at each time-step [71]. The study developed Artificial Neural Networks (ANNs), where gradient descent was responsible for minimizing the square sum of the difference between real outputs and network outputs. Other applications of such a method include heat load forecasting and efficiency maximization. In case of heat load estimation, this technique was used for the minimization of mean absolute percentage error associated with weather forecasting [72]. While the other study considered its functionality to solve the optimization problem and maximize the efficiency of a FC-CHP [73]. In 1996, one of the first attempts to employ the gradient descent method in a CHP system was taken by M. I. Henwood et al., in order to solve the economic dispatch problem [74]. The suggested algorithm was comprised of two layers, in which the outer layer was responsible to solve power dispatch using the Lagrangian technique; while in each iteration, the inner layer would solve the heat dispatch problem considering the unit heat capacity passed by the outer layer.

Constrained Optimization Techniques
Over the past decade, the field of constrained optimization has evolved rapidly, where today, it covers a vast territory of techniques and applications. Using the evolutionary principles for practical problem solving, and development of evolutionary algorithms made the biggest contributions to this trend. The term evolutionary computing is used to describe the field of analyses that involves heuristic search principles inspired by natural evolution [77]. Evolutionary computing evolved the whole concept of classical optimization, and it has been receiving increased interest by offering a variety of practical advantages to the energy models; facing problems concerning optimization, model identification, and control in power systems. The advantages being simplicity of the method, robust response to dynamic changes, flexibility of development, and broad applicability. In this context, a systematic literature review was conducted to gather available research articles employing evolutionary computing along with integer programming in IESs. The systematic literature review allows one to identify, select, and overview the preliminary research works thoroughly. In this regard, preliminary studies were identified and selected based on the use of primary keywords, concerning the implemented optimization technique, suggested application, and the proposed energy model. The studies are classified based on the implemented algorithm into four groups, including Genetic Algorithm (GA), Particle Swarm Optimization (PSO), Simulated Annealing (SA), and Mixed-Integer Linear Programming (MILP).
In this section, the four techniques are discussed precisely, and their application in constrained optimization of CHP and CCHP systems are explained. Further, statistical reports concerning the primary application of each method are provided as well. Improving the structure and performance of energy systems is associated with various restrictions imposed by technical limitations, circumstances, or organizations. Application of constrained optimization in energy systems include, but are not limited to thermal unit commitment, economic dispatch, optimal power flow, scheduling and planning, and improvement of thermodynamic performance. In general, such applications can be classified into four primary groups namely thermoeconomic optimization, thermodynamic optimization, economic optimization, economic dispatch, and scheduling. It is worth mentioning that environmental assessments are not considered in this classification since their solo evaluation (environmental-only assessment) is not contemplated in many studies. Often, environmental analysis and optimization are coupled with economic and/or thermodynamic assessments. Hence, in each category there might exist studies that considered environmental optimization along with other objectives. Similarly, other applications with nominal contributions, such as unitcommitment or load forecasting, are considered in the scheduling group. The systematic literature review illustrates that GA and MILP are the two methods with the most research contributions in the field of IES ( Figure 5). Further analyses on application of each method presents that, economic and thermoeconomic performance of CHP and CCHP systems are still of interest to the research community ( Figure 6). Thermoeconomic improvement of CHP and CCHP systems has been studied over the past decades, and with energy transition in action it will continue to be one of the main topics of research for the following years. Practical limitations such as economic dispatch, unit commitment, and scheduling in CHP and CCHP systems remain as complex problems, which require intensive mathematical formulation and computational programming. With improvements in evolutionary computing and nonlinear programming, such problems will be the research focus of future studies.

Genetic Algorithm Optimization
The Genetic Algorithm was first developed in 1989 by John Henry Holland at the University of Michigan [78]. This well-known method of optimization is also referred to as a stochastic search technique that is based on the mechanics of natural selection and genetics. This algorithm operates on a population of individuals to search several areas with possible solution(s), simultaneously and adaptively. In todays' power system, the application of this method is well recognized and it has been the most widely used technique not just for optimization, but also for system identification and control purposes as well [79]. The process of Genetic Algorithm comprises of four simple steps, including population generation, fitness function and selection, crossover, and mutation (Figure 7) [53]. In the first stage, an initial set of values that is guaranteed to include the final solution is randomly or heuristically created and named the initial population. The population specifies the upper and lower bound of search space and it contains the potential solutions to the given optimization problem. The solutions are known as chromosomes and shaped as fixed-length binary strings, comprising of 0 s and 1 s. In the second stage, the fitness function is applied to each potential solution to check how well the solution is responding to the problem. Then, the selection procedure is utilized to create an intermediate population with better suited chromosomes. In the following steps, crossover and mutation procedures are implemented in intermediate population in order to exchange some portion of the strings, and change the string locally so as to create a better string with a small probability. The entire process is carried out until the termination criteria are satisfied. Another interesting aspect of GA is its performance in handling constraints. This method is capable of processing constrains in two different ways, including embedding and penalty functions. While in most efficient cases the constraints are taken into consideration by embedding them in the coding of the chromosomes, in other optimization models each limitation is defined in terms of a penalty function to avoid ideal performance [80]. In general, the GA method of optimization presents various advantages over classical techniques. The primary advantages and disadvantages of GA technique can be described as follows:  Although this method provides undeniable benefits over the classical approach, but the first-generation models of this technique are incapable of solving complicated problems due to a large amount of computational time and process [81]. Employment of first-generation GAs in large-scale power systems has shown a dramatic increase in execution time with low quality solutions (random initialization) [82]. A selection of preliminary studies, based on research impacts and the optimization model, is presented in Table 1. The studies are classified based on the objective functions and type of constraints, in order to realize the current research focus and the primary application of GA in IES optimization. Evidently, majority of studies considered economic and Primary Energy Saving (PES) assessment of CHP or CCHP systems. Almost all studies attempted to minimize the Annual Total Cost Saving (ATCS) of the energy systems, while taking energy balance and design parameters constraints into consideration as well. Broad applicability of GA allows such a technique to be utilized for thermodynamic, economic, and environmental assessments, and coupling of these three perspectives. In thermodynamic optimization, often the objective(s) is to improve the thermodynamic performance with respect to the overall and exergy efficiency, hence, the optimization is limited to thermodynamic improvement. In this regard, the constraints are mainly evolved around components' limitation, energy balance, and design parameters. For instance, temperatures, flowrates, component's efficiency, and power-tohear ratio are the key constraints in every thermodynamic optimization model [93]. On the other hand, economic constraints are the major restriction of the economic optimization. In such optimization models, the objective functions are mainly evolved around economic performance, however, energy balance constraints may still exist in the optimization model. CHP or CCHP systems without heat and cooling balance constraints are just another power generation asset. To this end, economic constraints such as financial availability, maximum payback period, and legal constraints (related to governmental policies) are imposed along with energy balance constraints [94]. Economic dispatch can also be interpreted as economic optimization, in which the main objective is to minimize the production cost by properly scheduling the operation of the system. In the economic dispatch problem, energy balance, including power, heat, and cooling, and capacity limit, including power-only, heat-only, and CHP unit should be considered in the optimization model [95]. In economic dispatch, another important constraint is the Prohibited Operating Zone (POZ) in which the generating unit is restricted to its physical limitation. It is a state between maximum and minimum generation limits, in which the generating unit is not available due to downtime or fault in the system [96]. Several studies on CHP and CCHP economic dispatch considered energy balance as the primary constraint in the optimization model [97,98], while POZ was considered in only a handful of studies [95,99]. In the scheduling application, in addition to the energy balance, ramp rate and minimum ON/OFF time can also restrict the optimization model as well [100]. The scheduling model can be extended further by considering transmission constraints such as the three-phase electric network, resource availability constraints such as the natural gas network, and grid constraints such as the minimum load in the optimization model [89,101,102]. Environmental constraints can be presented and coupled with other thermodynamic and economic constraints. Common environmental constraints include legal CO2 emissions limitation, fuel consumption (depending on the perspective it can be interpreted as thermodynamic or economic constraints as well), and fuel quality [103,104]. In general, applicability of GA in optimization of CHP and CCHP is vast, and it can include a variety of constraints in the optimization model. Although, GA does not scale with complexity (large data set), therefore, it may not be the best solution for complex problems such as Combined Heat and Power Economic Dispatch (CHPED). In this regard, employment of such technique for scheduling application is limited, unless the CHPED is broken down into several simpler representation. In this case, GA can be employed as a suboptimizer for solving a certain objective function. Table 1. Selected preliminary studies on energy system optimization using the GA method.

References
Description Objective Function(s) Type of Constraint(s) [105] Analyzing the influence of CCHP on residential and office buildings. Optimization the proposed CCHP system.

PES ATCR CDER
Energy balance Tariff policy [106] Optimization of a CCHP system coupled with ground source heat pumps and PV/T under different load ratios.

APESR ATCSR ACDERR
Energy balance PV/T installation area [107] Optimization and analyses of a CCHP system comprised of supercritical CO2 recompression cycle and gas turbine waste heat.

TPUC TAPR
Equality and equality [108] Optimization of a hybrid CHP energy system integrated of solar, wind, and FC technologies for residential demands. Optimization of three different types of hybrid CHP energy systems considering cumulative energy demand.

PEC CED
Installation costs Architectural constraints [110] Optimizing the operation of a HT-PEMFC-CHP for residential sector under steady-state conditions. Net power output Net heat output Electrical efficiency Thermal efficiency Design parameters such as steam-to-carbon ratio and anodic stoichiometric ratio [111] Optimizing the operation of a micro-CCHP system considering energy coupling analysis PESR Energy cost Equality and inequality [112] Developing, evaluating, and optimizing a CAES-CCHP system Round trip efficiency ATCS Unit capacity and volume [113] Performing energy analysis on a CAES-CCHP system integrated with solar energy and optimizing its operation.

Exergy efficiency TIC
Design parameters such as thermal oil temperature and air compressor pressure ratio [114] Assessing and optimizing the operation of a CCHP system integrated with solar and geothermal energies from energy, economic, and environmental perspectives.
PESR ACSR CDERR Equipment capacity Energy balance

Particle Swarm Optimization
In early 1990, the research community began to realize the natural creatures' behavior and their applicability in different sciences including mathematics. This soon started to attract attention in the world of computing and algorithm. About the same time as M. Dorigo et al. [115] developed Ant Colony Optimization (ACO) based on social insects [116], J. Kennedy and R. Eberhart were formulating the Particle Swarm Optimization (PSO) algorithm based on the analogy of swarms of birds as they search for food [117]. This is a stochastic global optimization method, which is based on the exchange of information between individuals (particles), unlike creating a new population in GA. As the researchers examined the behavior of decision making in humans, they found out that all decisions are made based on their own experience and experience of others. The same concept is utilized in PSO formulation as well, where each particle in the population (swarm) is treated as a point in the 2-dimensional space with a certain velocity (Figure 8). The particle modifies its position based on its own experience (local best) and the experience of others around it (global best). This phenomenon allows particles to have exploratory and exploiting behavior so as to search and reach the optimal solution in the space through modification of particle's position and velocity. Similar to the GA method, each particle is considered as a possible solution to the optimization problem. Further, the fitness function is applied to each particle in order to identify four main characteristics of each particle including current position (x, y), current velocity (vx, vy), current distance between its position and its local best, and current distance between its position and the global best [118]. At every iteration, the velocity and position of each particle is modified and the performance is evaluated with respect to the fitness function and inertia weight, which indicates the impact of previous velocity and position on the current one. The process is carried out until satisfying the termination criterion by achieving the minimum distance between the current position and target or reaching the maximum number of iteration specified [53]. The PSO method is characterized with respect to three main factors including inertial weight, memory, and cooperation [53]. The primary advantages and disadvantages of PSO technique can be described as follows: Lower maturity level and commercialization compared to GA. Amongst the influencing parameters on the performance of PSO algorithm, inertia weight plays a significant role, reminiscent of the temperature parameter in Simulated Annealing (SA). This value determines the ability of PSO in both global and local search levels, and should be adjusted to balance both local and global exploration ability of PSO. Inertia weight value identifies the suitability of a suggested PSO for global or local optimizations. While the global search requires a higher value of inertia weight, the local search performance is negatively affected with higher values. To this end, most studies tend to develop adaptive PSO, in which inertia weight can be adjusted in accordance with the search space. Often, its value is higher at the beginning of the search, whereas, it shrinks as the solution gets closer to its optimal value [119]. Comparing to the GA technique, both methods start the process by initializing a population of candidates that contains all the possible solutions to the optimization problem. While GA is discrete in nature (binary variables), PSO can handle continuous problems. Although, both methods can be subjected to modification so as to handle discrete and continuous problems. The PSO search technique is formulated in a way so that all the members of the population follow the particle with the global best position value, whereas, GA requires information sharing between the chromosomes. In GA, the basic operators, including selection and crossover, are responsible for convergence rate modification, whereas, velocity equations are the ones for PSO. Since PSO is defined in 2-dimensional space, this technique is more suited for problems with less than three parameters, while in contrast, more than three parameters are required for GA. In this regard, PSO takes advantage of simpler structure and implementation, which also provides less computational time and less allocated memories [120,121]. Although, this is only guaranteed if the inertia weight is selected properly (adaptive PSO) [12]. First-generation PSO algorithms, similar to other heuristic optimization techniques, suffer from lack of mathematical foundation and they are likely to fall into local optimum in high-dimensional spaces. Moreover, PSO algorithms are faced with limitations in real-time economic dispatch problems due to their high computational time compared with other gradient-based techniques. Nonetheless, PSO is considered as one of the most powerful tools for optimizing mixed-integer non-linear problems [122]. Its simplicity and robustness make this method of optimization as one of the most utilized techniques in todays' power systems. Considering its undeniable advantages in energy systems optimization, researchers started to modify and improve this method in order to develop an algorithm based on PSO with better suitability for certain applications or variables. Today, there exist at least seven different variations of this technique, namely:  Discrete PSO;  PSO for MINLPs (optimal scheduling of CHP and manufacturing [123]);  Constriction Factor Approach (CFA);  Hybrid PSO (solving the economic dispatch problem of a CHP system [124]);  L-Best Model;  Adaptive PSO (solving the economic dispatch problem of a CHP system [63]);  Evolutionary PSO (developing and optimizing a hybrid-CHP [125]).
The statistical report of PSO developments in co-and tri-generation systems indicates strong interest in this method for economic optimization ( Figure 6). This can be due to the fact that the economic and dispatch problems have complex and nonlinear characteristics with various constraints, which makes PSO an interesting option for solving such problems [126]. Objective functions and operating constraints are relatively similar to the ones that are discussed in Section 3.1, due to the fact that such factors are objective-based rather than method-based. In accordance to the summary provided in Table 2, operating cost and total cost of production are the primary objective functions in the optimization models. The majority of studies considered equality and inequality constraints along with design parameter limitations such as temperature and expansion ratio in the turbine [127]. In economic dispatch problems, constraints could include ramp rate and minimum ON/OFF period, such as in Refs. [128,129]. Economic constraints can be presented for a single component or an input to the system, such as the cost of working fluids in Organic Rankine Cycles (ORCs) in Ref. [130], and maintenance cost of reciprocating engine in Ref. [131]. Similarly, the equality and inequality constraints can be introduced for a single component, such as the pumped storage unit in Ref. [132], to limit the physical size and operating conditions to the realistic boundaries. With regards to the CHPED problem, various studies incorporated POZ constraints into the optimization model [63,133,134]. Concerning thermoeconomic optimization, improvements in mathematical optimization attracted the attention towards Linear and Non-linear Programming (NLP); due to the fact that such optimization techniques can yield the lowest percentage of gap, compared to heuristic methods. To this end, contributions of PSO to thermodynamic optimization have been limited, since mathematical methods such as MILP can obtain the most accurate results. Although the implementation of PSO has not been investigated in many studies, compared to GA, the current developments clearly shows that PSO is still under investigations and improvements, and one can hope for its vast implementation in CHP or CCHP systems in the near future. Table 2. Selected preliminary studies on energy system optimization using the Particle Swarm Optimization (PSO) method.

References Description
Objective Function(s) Type of Constraint(s) [135] Optimization of a hybrid solar-hydrogen based FC-CHP energy system with two heuristic methods using adaptive inertia weight and constriction factor Overall operating cost Overall income Design parameters such as rated capacity of FC Ramp rate [136] Developing an off-design optimization model for a CCHP-GSHP from thermos-economic and environmental point of views  [142] Considering the environmental emissions cost in allocation of renewable and nonrenewable CHP systems.

Overall operating cost
Overall income Design parameters such as Header pressure tolerance

Simulated Annealing Optimization
Simulated Annealing (SA) is known as the most flexible and preferred heuristic method for solving complex combinatorial problems. The approach was inspired by an annealing procedure of metal working, which defines the solidification and formation of crystals by slow cooling of substances to the minimum energy state [143]. The most important advantage of SA over classical or other heuristic methods is the capability of handling large problems regardless of variables' conditions. SA approach was formulated when the relation between the physical state of the matter and the solution space of an optimization problem was realized. The principle of this analogy states that an optimal solution is associated with a perfect crystal, while a local optimum indicates a crystal with defects [118]. The SA method is based on a point-by-point search technique, which utilizes a control parameter that has to be properly determined for achieving efficient results. The simplest type of this algorithm compares iteratively the output of the objective function, considering the current and neighboring point in the domain. If the neighboring point generates a better solution, then it is considered as the base solution for the next iteration, otherwise, the procedure is terminated without searching the wider domain. This helps the algorithm to avoid convergence to a local minima or maxima. Instead of widening the search domain, SA takes advantage of two iterative loops, namely the controlled cooling procedure and Metropolis criterion. In the optimization process, the algorithm attempts to minimize the value of loss function by capturing the procedure of controlled cooling mathematically. Although, in the face of the local optimum, the value of the loss function might experience a provisional increase since it gathers necessary information to reach the global minimum, which is accounted for in the SA model. Similar to the adaptive inertia weight of PSO algorithm, a control parameter (known as temperature) is introduced based on the concept of Boltzmann probability distribution, which is varied over time using nested loops. At the beginning of the search, the temperature is set at a high value allowing larger deterioration in the cost function ( Figure 9). While the temperature is incremented and the search sector is narrowed, the algorithm becomes greedier and only accepts smaller deteriorations. At every iteration, the Metropolis criterion allows the process to be executed randomly in order to extra search the neighborhood so as to avoid the local extreme points. The new state should satisfy the Metropolis criterion, which has the same form as Boltzmann's Equation [144]. Through its nested loop, it adopts the local minimum at each inner loop, and correspondingly, it helps to eliminate such points from the search space. To this end, this technique is amongst the most preferred heuristic approaches that provides practical randomness into the search. The control strategy of SA is based on the principle of Cooling Scheduling in the crystallization operation. It is used from the starting point until the convergence is reached. The primary advantages and disadvantages of SA technique can be described as follow: Large data sets can be computationally too expensive. Surprisingly, not many studies contributed to employment of SA in co-and trigeneration systems, from which mostly elaborating on the economic aspects of CHP and CCHP configurations ( Figure 6). This can be due to the fact that such a method is overskilled for optimization problems related to thermodynamic aspects of IESs, and there exist simpler methods such as PSO and MILP to solve for such problems. The SA technique is mainly employed for combinatorial optimization problems, which makes this method an interesting choice for scheduling, unit commitment, and economic dispatch problems in power systems and IESs [145]. The majority of contributions to CHP and CCHP systems are also in the forms of scheduling and economic improvement for microgrid and smart-grid applications. In Ref. [146], an improved version of SA was employed in order to find the optimal energy management strategy in a microgrid network with renewables integration and possibility of distributed generation. The objective of optimal energy management strategy in microgrid networks is maximum financial gain, while prioritizing the distributed generation using renewable energies. In a similar study, SA was utilized for day-ahead resource scheduling in a smart-grid network [147]. The objective function was minimization of the Virtual Power Player (VPP), while taking into account various constraints such as energy balance, AC voltage balance, voltage transformer capacity, and storage capacity. A general overview of preliminary studies is presented in Table 3, which indicates the importance of the SA method in the economic analysis and optimization of power systems. Majority of studies considered fuel cost of the system as the objective function to be optimized. Given the lack of available contributions to IESs, the presented preliminary studies are not bounded to CHP and CCHP systems. In fact, the application of this optimization technique was considered in a larger scale to include power systems in general. Realizing the functionality of this technique might help the researchers in the field to take interest in this method of optimization and investigate its feasibility for energy integration, particularly CHP and CCHP systems. Considering the advantages of SA in economic improvement, its implementation as a suboptimizer in a framework along with PSO, GA, or even Linear Programming (LP) can be a great approach towards developing hybrid optimization algorithms and solving the CHPED problem. Table 3. Selected preliminary studies on energy system optimization using the SA method.

References
Description Objective Function(s) Type of Constraint(s) [148] Developing and optimizing a solar/wind hybrid energy system subjected to constrained TAC CRF Equality and inequality Number of PV panels Number of wind turbines [149] Developing a fluid-selection and optimization algorithm using SA for an ORC system Heat exchanger area to net power output Net recovery efficiency Design parameters such as condensing pressure [150] Developing an algorithm based on PSO and SA for predicting bidding strategy considering lack of information about market mechanism ISO's objective function Equality constraints Transmission line constraints Generation capacity constraints [151] Developing and optimization a hybrid energy system based on solar and wind energies LPSP PBP Number of PV panels Number of Battery units [152] Solving the short-term UC problem of a CHP system using SA integrated with dynamic economic dispatch method Fuel cost Start-up cost Power balance constrain Ramp rate [153] Developing a MG planning model considering optimal sizing of the hybrid CHP system AC Design parameters such as installed capacities Operating constraints Transmission constraints [154] Formulating a unique algorithm based on SA and EA method for solving the economic load dispatch problem for a power system Fuel cost function Inequality constraints Power balance constraint [155] Solving dynamic economic emission dispatch by proposing a unique algorithm based on SA principle of operation Fuel cost Emission rate Equality and inequality [156] Solving the combined power and economic dispatch of a power system considering price penalty factor Fuel cost function Equality and inequality [157] Developing an algorithm based on ABC and SA for solving the combined economic and emission dispatch considering constrained and non-linear multi-objective function Fuel cost Total emission Power balance constraint

Mixed-Integer Linear Programming Optimization
Over the past decade, Mixed-Integer Linear Programming (MILP) formulation has gained considerable prominence in the field of energy generation and optimization. In many optimization problems, integers are the only appropriate values for system variables; this type of problem formulation is known as Integer Programming. Such a method of problem formulation is one of the most common and powerful mathematical modeling techniques for a variety of engineering systems with different optimization problems including, but not limited to: unit commitment, economic dispatch, and thermoeconomic assessment. Objective functions and decision variables can take up linear and non-linear forms, which can affect the use of Integer Programming in a system. Undoubtedly, linear problems are simpler to understand and require less computational efforts to be solved; therefore, developers prefer LP over the NLP form. In recent time, there has been seminal contributions in the field of Non-Linear Optimization. The comparative analyses between LP and NLP have shown the superiority of NLP performance with regards to the accuracy and reliability of the optimal solution [158]. Although, the NLP technique is associated with major downsides, which make it unreliable for certain applications in energy systems. Some of the disadvantages are as follows:  Converging to a local optimum is very likely;  Initialization affects the final solution considerably;  Difficult to satisfying equality constraints;  Complexity of formulation and modeling;  Combinatorial explosion of integer variables.
To this end, energy systems are mostly developed in linear forms, and non-linear functions are subjected to linearization methods for further simplification. Considering any engineering structure, there exist variables that do not require integer values, or the formulation is not constrained to use discrete variables; hence, the model can be developed with both integer and non-integer values. This method of formulation is referred as Mixed-Integer Linear Programming, which is classified as a discrete optimization problem [159]. Despite the existence of the MILP approach over the decades, this method has been significantly developed and is now one of the simplest, widely used, and robust optimization tools in energy and power engineering. The structure of this programming method is based on linear functions, however, there is no restrictions over the use of fractional values in the programming. In this context, objective functions or constraints are able to take discrete or continuous forms in MILP models. The problems formulated with MILP are solved using the branch and bound algorithm, in which the problem is branched into nodes and then solved as relaxed Linear Programming problems [158]. The branching procedure is continued until the required tolerance between the integer solution (upper-bound) and Linear Programming solution (lower-bound) is reached. With current developments and available methods of optimization, the MILP approach is the best possible option to consider for the formulation of power systems in comparison with non-linear programming. It is clearly a common language to uniquely describe a problem in strict mathematical terms. Although, this does not indicate that such a technique is not associated with any drawbacks. The primary advantages and disadvantages of MILP technique are described below. It is worth mentioning that such downsides can be overcome by properly formulating the optimization problem. There are a variety of methods for solving MILP problems, namely: the Branch-and-Bound Method (economic optimization [94]), Cutting Plane Method (load forecasting and dispatch optimization [160]), and Branch-and-Cut Method (life-cycle optimization [161] Expectedly, this technique has been the second most implemented method for IES optimization. Given the mathematical nature of this technique, the applicability of MILP in economic dispatch optimization is vast. Readers are advised to refer to Ref. [162] for further information on practical applications of MILP, and linearization of nonlinear problems. Commercialization of MILP through numerous platforms such as MATLAB and CPLEX has provided supporting managerial decisions, and increased interest in employment of this technique for complex nonlinear problems. In CHP and CCHP systems, vast majority of employments contributed to economic improvements, particularly hybrid renewable-based systems, so as to encourage the private and public sectors for further consideration. With current developments in the IES field, this method of optimization along with MINLP will be the most utilized techniques in the near future. A summary of selected studies is presented in Table 4. It can be illustrated that the total cost, operating and maintenance cost, and annual profit are the primary objectives functions in most of the studies. The optimization models are mainly subjected to a variety of constraints including, but not limited to, energy balance, operating conditions, and design parameters. There have been nominal contributions on micro-CHP systems for residential applications in order to utilize MILP for optimal capacity and the dispatch profile formulation [163]. In some cases, environmental objectives functions, such as minimization of annual GHG emissions, were considered along with economic ones [164]. Multiobjective economic dispatch optimization of CHP and CCHP systems is one of the topics with the highest contribution. Many studies considered several objective functions such as revenue, generation cost, storage cost, and cost of firm and non-firm power exchange, in which the optimization models were restricted to various constraints such as energy balance, network constraints, charge and discharge rate of storage units, voltage balance, and grid security [165,166]. A handful of studies formulated the economic dispatch and scheduling problems into a two-stage (or levels) stochastic model, in which every stage can be responsible for a certain objective. For instance in Ref. [167], the commitment of the electricity and heat generation units with respect to the power and heat demand would be determined in the first stage, while the operating and production costs would be determined in the second stage appropriately. A similar approach was taken in different studies with diverse energy system models to coordinate an optimal day-ahead planning, considering different constraints [168,169]. Since MILP is not a search-based technique, it does not follow a precise procedure. To this end, an example of MILP implementation for scheduling application in smart buildings utilizing CHP units is presented in Figure 10 Table 4. Selected preliminary studies on energy system optimization using the MILP method.

References
Description Objective Function(s) Type of Constraint(s) [171] Proposing a technology selection and operation optimization model for distributed energy systems considering cost and emission factors TC Total emission Total cost of system and emission Technology selection constraint [172] Suggesting an optimal dispatch strategy for CCHP system integrated with wind power considering operating constraints Developing a dispatch model for coordinated operation of an Energy Hub consisted of CHP unit considering environmental and economic aspects TC Equality constraints Power system constraints NG system constraints [178] Developing a novel operating scheduling for a CHP-based IES considering thermal inertia of buildings and renewables uncertainties Generation cost Energy balance constraints Thermal power unit Components constraints [179] Proposing a self-healing in distributed generation with gas and electrical network, considering various types of technologies such as CHP

TOC
Gas and electrical network constraints [180] Developing a thermoeconomic optimization model for energy, reserve, and reliability provided by a MG. Considering Transactive Energy approach for co-optimization purposes TC Energy balance constraints Reserve and reliability services

Conclusions
Design and development of energy systems can be classified into various steps, including, but not limited to, planning, model development, and forecasting, amongst which the optimization is one of the most important steps. Optimizing the performance of an energy system is improving the behavior and characteristics under steady-state and extreme conditions in order to make the most effective use of resources while achieving the primary goals. The optimization techniques have different natures and operating principles. While mathematical methods make use of derivatives to find the precise optimum solution, heuristic techniques search for the best possible solution amongst a large population of candidates at each iteration in order to reach the full convergence.
CHP and CCHP systems are two of the most well-known and fully developed energy integration concept. Such systems are capable of simultaneous production of power, heat, and/or cooling in the priority or form that is required. This type of integrated energy system is associated with many advantages such as fuel diversification, lower environmental impacts, and improved economics. Over the past decade, CHP and CCHP systems have attracted interests particularly for distributed generation application in the industrial sector. Although, the application of CHP and CCHP systems can include centralized generation if thermoeconomic and dispatch problems are to be solved. In this context, various optimization techniques were reviewed in this study that could improve the thermoeconomic performance of such systems. Depending on the availability of constraints in an optimization model, the methods were classified into constrained and unconstrained techniques. It was discussed that the availability of constraints would significantly add to the complexity of CHP and CCHP systems. To this end, employing evolutionary algorithms could help to optimize the performance using heuristic techniques such as PSO or GA, or taking advantage of linear mathematical programming with MILP. On the other hand, there exist various conventional unconstrained techniques that are still robust and reliable to use for optimization of energy systems. These techniques are associated with less complexity in comparison to the evolutionary algorithms. In terms of application, the MILP method was mainly utilized for economic optimization while GA and PSO were employed for a wide range of applications including thermodynamic, economic, and thermoeconomic improvements.

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