Optimization, Stability, and Entropy in Endoreversible Heat Engines

The stability of endoreversible heat engines has been extensively studied in the literature. In this paper, an alternative dynamic equations system was obtained by using restitution forces that bring the system back to the stationary state. The departing point is the assumption that the system has a stationary fixed point, along with a Taylor expansion in the first order of the input/output heat fluxes, without further specifications regarding the properties of the working fluid or the heat device specifications. Specific cases of the Newton and the phenomenological heat transfer laws in a Carnot-like heat engine model were analyzed. It was shown that the evolution of the trajectories toward the stationary state have relevant consequences on the performance of the system. A major role was played by the symmetries/asymmetries of the conductance ratio σhc of the heat transfer law associated with the input/output heat exchanges. Accordingly, three main behaviors were observed: (1) For small σhc values, the thermodynamic trajectories evolved near the endoreversible limit, improving the efficiency and power output values with a decrease in entropy generation; (2) for large σhc values, the thermodynamic trajectories evolved either near the Pareto front or near the endoreversible limit, and in both cases, they improved the efficiency and power values with a decrease in entropy generation; (3) for the symmetric case (σhc=1), the trajectories evolved either with increasing entropy generation tending toward the Pareto front or with a decrease in entropy generation tending toward the endoreversible limit. Moreover, it was shown that the total entropy generation can define a time scale for both the operation cycle time and the relaxation characteristic time.


Introduction
The optimization of energy converters has never been as relevant as it is now. Energy production requirements, efficient use of heat sources, and heat waste reduction are continually pushing technological edges. Linked to this degree of specialization are the control and stability of operation regimes yielding a desirable stable production, despite the possible variations of external or internal conditions. In many cases, this requires the fine-tuning of control parameters with intrinsic energetic costs. In this regard, there are some hints as to the possibility of seizing the stability of heat engine operation regimes to enhance their performance by relaxing the control over the operation parameters [1][2][3]. Studies regarding the weakly dissipative limit [4] thus far have pointed out the special role played by the endoreversible model [5,6]. These studies consider that after the system experiences a perturbation in its operation variables, the trajectories that lead the system back to its steady state tend to evolve toward the endoreversible limit. This limit has been proven to be associated with thermodynamic states, with the best compromise between maximum efficiency, maximum power

A Quick Look at the Endoreversible Model
The most basic endoreversible model consists on a baseline Carnot cycle whose working fluid operates with effective temperatures T hw and T cw (isotherms of the Carnot cycle), irreversibly connected to two external reservoirs at temperature T h and T c (T h > T hw > T cw > T c ), as depicted in Figure 1. The endoreversible scheme requires knowledge of the heat transfer laws to model the heat fluxes between the external reservoirs and the working fluid. In a somewhat general case of heat transfer laws, Q h and Q c are expressed as: where k = 0 is the exponent of the heat transfer law (k = 1 refers to the known Newton law, k = −1 is known as the phenomenological law, k = 4 is the Stefan-Boltzmann law, and so on); σ c and σ h are the thermal conductances with units of J·K −k /s, here considered as constant since only small fluctuations around the operation regime are addressed. The function sgn (k) is defined as: which is introduced for consistency with the convention that the heat flux entering the working fluid is positive and that the outgoing heat is negative; t h and t c are the times at which the working fluid is in contact with the external heat reservoirs T h and T c , respectively. In Equation (1), the σ hc ≡ σ h /σ c ratio is introduced, which is an important ingredient to obtain the upper and lower bounds of the efficiency of a given operation regime. The endoreversible hypothesis states that the entropy generation in the inner part of the engine is zero. Mathematically, this means that: and from this constraint, the ratio of the operation time is given as: where τ ≡ T c /T h . With the endoreversible hypothesis constraining t c , in Equations (3) and (4), the efficiency η, power output P, and total entropy generation S become functions of {T h , T c , T hw , T cw , σ hc , t h , σ c }. By using Equations (1) and (2), the explicit equations of these key thermodynamic magnitudes are the following: As can be seen in the definition of the power output, instantaneous adiabatic processes were considered. In order to deal with non-instantaneous adiabatics, a more detailed analysis of the working fluid and the geometry of the device needs to be made (for example, see [36]). In some works, these features have been analyzed by showing their influence in the power output and efficiency, but under certain circumstances, such as the large compression ratios in a piston, these times produce small effects on η and P, remaining compatible with the quasistatic nature of the internal process.
Notice that only S is proportional to t h . As is shown later in Section 4 (see Equation (15)), the relaxation times are also proportional to this partial time; thus, a relationship between relaxation time and total entropy generation can be found. By fixing the value of S, a time scale for stability can be established. Besides this time scale (establishing the speed of the restitution dynamics), it is possible to analyze the energetic consequences of stability independently of t h if one works with entropy production per cycle time, defined asṠ ≡ S/(t c + t h ) = S/t h (1 + t c /t h ), taking advantage of the constraint on the total operation time through Equation (4).
The maximum power (MP) regime is specified by the temperatures T * cw resulting from the constraint ∂P/∂T cw = 0 and T * hw under the condition ∂P/∂T hw = 0. The resulting values are given by: From the optimal T * cw values, parametrization of P, η, andṠ can be made (see dashed curve in Figure 2). The obtained parametric curve exhibits a parabolic-like behavior, where the values of η vary from 0 to η C ≡ 1 − τ, the Carnot efficiency (both limits corresponding to a zero power output), and with a single maximum in P. This is what is known in the literature as an endoreversible behavior, a kind of signature of the model. Some relevant features regarding the efficiency at maximum power from these endoreversible models with k = 1 and k = −1 are pointed out:

•
The Newton case (k = 1) gives an efficiency at maximum power, η MP = η CAN ≡ 1 − √ τ, the well-known Curzon-Ahlborn-Novikov (CAN) efficiency, which does not depend on the σ hc ratio and appears in a large variety of contexts linked to the maximization of power output, work, and kinetic energy [37].

•
The k = −1 case, frequently called the phenomenological law, referring to the natural results arising in the linear irreversible thermodynamics framework. It allows to obtain the same limits of efficiency as in the low-dissipation model, where the self-optimization property has been studied.

The Relevant Region for Optimization: The Pareto Front
Before exploring stability dynamics, it is convenient to introduce a multiobjective optimization of the model. When looking for the best compromise among a variety of objective functions, the result is the so-called Pareto front (in the space of energetic functions) and the corresponding Pareto optimal set (in the space of operation variables).
To this end, a sorting algorithm is used by applying the concept of dominance [38]: A vector v = (v 1 , . . . , v n ) dominates another one w = (w 1 , . . . , w n ) if, and only if, v i ≥ w i ∀ i ∈ {1, . . . , n} (if one is looking for a maximum, ≤ for a minimum) and there is at least one j, such that v j > w j , if one is interested in those vectors that are not dominated by any other. In other words, vectors that have the best value at least in one objective function. In this case, such a vector is formed by η, P, andṠ. The algorithm introduced here is a modification of the one introduced in [1-3], as follows: 1.
In the phase space (T hw , T cw ), the region of physical relevance is defined (T h ≥ T hw ≥ T cw ≥ T c ).

2.
A random set of points in the phase space is obtained and the thermodynamic functions are evaluated (energetic space).

3.
A set of non-dominated points in the energetic space is obtained, giving a provisional Pareto front.

4.
From the corresponding Pareto optimal set (phase space), a convex region is defined and extended in order to cover a larger region for searching new points in the Pareto front. Details on the definition of the extended region are given below.

5.
From the new region, a new set of random points is proposed and a new set of non-dominated points in the energetic space is obtained.
In the present analysis, and in order to ensure convergence in the results, Kullback-Leibler (KL) divergence was introduced as a measure of the relative entropy D KL [39]. D KL was calculated between the probability distribution of the efficiencies of the Pareto front in the ith and the i − 1th iterations. As the probability distribution of η converges to the true Pareto front distribution, the entropy of the distribution converges as well. In such a case, the relative entropy D KL tends toward zero. The radius to extend the search region in the phase space decreases with the D KL value. When this relative entropy is very small, there is no information gain in iterating more times; then, the search for new points in the Pareto optimal set stops. As mentioned before, D KL provides a measure of statistical convergence by indicating how distant two distributions are. If D KL = 0, then the information stemming from both distributions is the same. This is a relevant issue to demonstrate that the obtained trend is not due to the lack of additional iterations.
To obtain D KL , the range of possible values of η, each iteration was divided into √ N (rounded to the upper next integer [40], with N being the number of random points added to the search in each iteration) equal intervals (or bins). In this way, the same partition was used to compute the discrete probability distributions, ρ k (for the k iteration). D KL was calculated by comparing ρ k−1 with ρ k . D KL,k is given by: allowing to determine how much information is gained by narrowing the search. In Figure 2, the Pareto front is depicted (green points), along with the endoreversible curve for η, P, andṠ (dashed line), obtained from the parametric elimination of T hw (T cw is constant for the fixed values of the σ hc , σ c , τ, and T h parameters; see Equations (7) and (9)). Notice that even when the endoreversible curve depends on a first constraint (T cw satisfies ∂P/∂T cw = 0), the Pareto front has nothing to do with. Nonetheless, the Pareto front lies over the endoreversible curve, covering the regions from maximum power to maximum efficiency and minimum entropy production. This is consistent with the literature, where the mentioned part of the curve has been denoted as the relevant region for optimization. For optimization purposes, additional figures of merit as compromise functions between these three (such as the Ecological function [41] or the Omega function [42]) will not provide further information to the Pareto front. After all of these considerations, analysis of stability dynamics can be made.

Stability Dynamics and Relaxation Times
One goal of this analysis was to recover the previous results obtained in the low-dissipation scheme (with no direct connection to working fluid particularities) and to lay the groundwork for a more direct connection with the properties of the working fluid and design parameters. The perturbations we were interested in involved only those of an external nature. The temperature of the external reservoirs remained constant, as well as the conductances of the heat transfer laws (lastly depending on temperature) for only small variations in the operation regime. To illustrate these points, consider a 2-D piston describing a Carnot cycle. The velocity of the piston will determine the effective temperatures of the isotherms of the particles inside the piston (see [32]). The external control (and its energetic cost) needs to remain constant as the velocity of the piston is not accounted for in the thermodynamic description of the gas inside of it. However, fluctuations in this velocity, as well as possible cyclic variability, will lead to variations in the effective temperatures T hw and T cw . Then, a compromise between the energetic cost of the control and the thermodynamic consequences due to stability is of interest.
A simple approach to tackle this problem is the following: Since the operation regime is entirely defined by T hw and T cw , the dynamics involving only these two variables are to be addressed. Fluctuations in these temperatures will affect the heat exchanged between the working fluid and the external reservoirs, Q h and Q c . From a Taylor expansion in the first order, Q h and Q c near the MP state can be written in matrix form as: where α, β, and γ are the elements of the Jacobian matrix evaluated in the MP state (denoted by * ). It is noted again that this analysis assumes that the MP state is a steady stationary state. Within the first order scheme [43], the simplest relationship for an autonomous system for the two variables T hw and T cw is: which is a good approximation near a stable point. The coefficients A and B determine the restitution strength (with units of s −1 ). The bigger they are, the faster the system will return to the steady state. By substituting Equation (12) into Equation (13), with the resulting dynamics: The magnitudes of the relaxation times t 1 = −1/λ 1 and t 2 = −1/λ 2 are obtained from the eigenvalues, λ 1,2 , of the square matrix in the right-hand side of Equation (14). For the cases k = 1 and k = −1, they are: Due to the units of the two matrices, a K/J unit factor should be accounted for in the final expression on the right-hand side of Equations (15), providing the correct units for the relaxation times (in s). As mentioned before, A and B provide the strength of the restitution dynamics, i.e., if they have large values, then the relaxation times are short. Notice that both the relaxation times and the total entropy generation are proportional to the contact time t h ; then, S can be used to define a characteristic time scale for relaxation. Since there are no reasons to provide a preferred relaxation in the heat exchange Q c or Q h , it can be assumed that t 1 = t 2 (this requirement can be tuned according to the specific conditions of a heat device at hand). Additionally, when looking for stability of a cyclic process, it is desirable that the stability is achieved in times shorter than those of the operation time; that is: From this constraint and Equation (4), the relaxation times are bounded by: where the equality corresponds to the case where t c + t h = t 1 + t 2 . From now on, equality is assumed. Then, the resulting dynamics (Equation (14)) for k = 1 are: and for k = −1, From these dynamic equations, it is possible to calculate a relaxation velocity, v = Ṫ 2 hw +Ṫ 2 cw (in K/s), which indicates how fast the system is evolving toward the steady state. In Figure 3, iso-velocity contours are depicted for the k = {1, −1} cases. Two trajectories are represented to provide an idea of how fast the system can return to the MP state. In both cases, the operation time is indicated, along with the time of evolution of each curve. As can be seen, 1 or 2 s are enough to drive the system close to the stationary state; meanwhile, the cycle time is approximately 300 s. Another feature in Figure 3 is that far from the steady state, the system evolves faster. Different marks at equal time intervals are placed over the trajectories, and in every case, the system travels a longer distance in the first interval. In both cases, the T cw direction is slower. It can be confirmed that the depicted trajectories evolve more quickly in the horizontal direction, while the final approach is mostly in the vertical direction. Notice also that the k = −1 case has the fastest dynamics. The values T h = 500 K, σ hc = 1, and τ = 0.4 were fixed. As in Figure 2, σ c and t h were chosen (for representation/comparison purposes) in such a way that the MP was 70 W and the entropy at MP was 70 J/K in both cases. By fixing S at MP conditions, a time scale for t h was established, and therefore, a scale for the relaxation time. The qualitative behavior for the other values of S was similar to the one presented here.

Thermodynamics of the Relaxation Trajectories
From the analysis of the evolution of trajectories toward the steady state, it is not obvious which energetic implications arise. To address this issue, the dynamic equations were numerically solved. Several trajectories with a starting point within a radius of 30 K (in the T hw -T cw space) from the MP state were analyzed. Then, these trajectories were mapped into the energetic space (η, P, andṠ). (4), (17) and (18), the conductance ratio σ hc plays a key role both in the dynamic equations and in the energetic properties of the system. In this regard, it is noted that the limits σ hc → {0, ∞} are not physical from a dynamic point of view. For both the k = {−1, 1} cases, the limit σ hc → 0 leads to null power output P, entropy S (Equation (5)), and also null restitution strengths (see Equation (18)), i.e., {P, S, A, B} → 0; thus, there is no power output and there are no restitution forces. In comparison, the limit σ hc → ∞ leads to {S, A, B} → ∞; thus, there is an infinite entropy generation and infinite restitution forces. For this reason, in the present analysis, the more realistic cases σ hc = {10 −6 , 1, 10 6 } were considered.

As it can be checked in Equations
In Figure 4, the four rows display the trajectories obtained from numerically solving the dynamics in Equation (18) for the k = 1 case. The first row displays the trajectories in the T hw and T cw space, while the next rows show, respectively, the mapping of the same trajectories over the η-P,Ṡ-P, andṠ-η spaces. The three columns are for the abovementioned representative conductance values: (a) σ hc = 10 −6 , (b) σ hc = 1, and (c) σ hc = 10 6 . In the first row, according to the initial state in each trajectory, a clock hour analogy was used; in this way, the purple line corresponds to 12 o'clock and red line to 3 o'clock. The Pareto optimal set is displayed (green points), together with some iso-velocity contours. The white curve indicates the position of the MP state as σ hc varies from 0 to ∞. All of these points produce the same efficiency at MP given by η * = η CAN . For the rest of the rows, the Pareto front and the endoreversible curve are presented as well. Note also that in the three configurations (columns), the total operation time is the same, but the t c /t h (and σ c ) ratio varies, as can be seen in the legend in each case. This figure reveals several key features: (a) For a small σ hc (left column), only trajectories between 9 and 3 o'clock are present in the T hc -T cw space (denoted in colors ranging from red to blue). Most of the observed trajectories (especially those with darker colors) evolve to the stable state near the endoreversible curve in such a way that the power and efficiency progressively increase their values, while entropy production decreases (see rows 2-4 in the first column). In other words, the relaxation process mainly drives the system toward a thermodynamic steady state, thereby enhancing the thermodynamic performance of the engine (η and P increase andṠ decreases).
(b) For σ hc = 1 (second column), the trajectories arriving at the stable point from all directions are clearly observed. Those in darker colors (from 9 to 3 o'clock) evolve, as in the above case, but closer to the endoreversible curve, showing an increase in power and efficiency and a decrease in entropy generation. In contrast, the trajectories denoted in colors ranging from orange to green (between 3 and 7 o'clock) evolve to the steady state near the Pareto front, with an increase in power output but a decrease in efficiency and higher entropy production. Note also in this configuration how the relaxation is well balanced between the trajectories approaching the endoreversible limit on the Pareto front side and to the other side. However, the curves near the Pareto front are longer, meaning that random perturbations tend to favor a locus directed toward the Pareto front.
(c) For a large σ hc (last column), the trajectories evolve directly toward the endoreversible curve first. In this part of the trajectory, η, P, andṠ are enhanced simultaneously. Later on, the system evolves to the steady state, either through the endoreversible limit or the Pareto front.
(d) Although the sizes of the perturbations in the T hw -T cw space are the same, the case with the small σ hc allows larger variations of power output, while the case where σ hc = 1 exhibits the smallest fluctuations in the η-Ṡ plane.
(e) The relaxation times (Equation (17)) are directly proportional to t h , and from the expressions of P, η, and S (see Equation (5)), only entropy generation depends on t h . Thus, the characteristic time scale of the relaxation is linked to the entropy scale of the system, a feature that connects the stability with the thermodynamics of the system. In this entropy-control point of view, the dynamics for the small values of σ hc in the left column apply to greater t h values and, as a consequence, when the contact time of the heat engine with the cold reservoir is small. On the contrary, the dynamics for the large values of σ hc in the right column apply to smaller t h values, i.e., when the contact time of the heat engine with the cold reservoir is large.  Figure 2, σ c and t h were chosen (for representation/comparison purposes) in such a way that the MP was 70 W and the entropy at MP was 70 J/K in both cases.
In Figure 5, the corresponding configurations for k = −1 are depicted. Notice that the behavior of the trajectories is almost the same as that in the k = 1 case, but here, the maximum power stationary state is linked to different efficiencies as σ hc increases from 0 to ∞. In particular, when σ hc → 0, the efficiency at maximum power is η * = 1−τ 2 = η c 2 ; when σ hc → 1, The dynamics for different σ hc are no longer linked to equal operation times t c + t h . Another difference is that the relaxation times are noticeably smaller and, therefore, the restitution forces are stronger. Additionally, as σ hc decreases (η * also decreases), the same temperature perturbations lead to the farthest starting points in the energetic planes. In this way, the case of tlarge σ hc (with a high efficiency) is the one with the smallest drops in η, P, andṠ (with shorter trajectories). This suggests that the system becomes more stable as the efficiency in the steady state increases.  Figure 5. Trajectories in the T hw -T cw space and mapping over the η-P,Ṡ-P, andṠ-η spaces for k = −1.
In column (a) the σ hc = 10 −6 case; (b) the σ hc = 1 case; (c) the σ hc = 10 6 case. The values T h = 500 K and τ = 0.4 were fixed. As in Figure 2, σ c and t h were chosen (for representation/comparison purposes) in such a way that the MP was 70 W and the entropy at MP was 70 J/K in both cases.

Concluding Remarks
It was shown that the Pareto front, which represents the best compromise among power output, efficiency, and entropy generation is related to endoreversible behavior, obtained herein by analyzing the Newton and phenomenological heat transfer laws in the context of Carnot-like models. These findings corroborate those obtained in the low-dissipation scheme.
A set of dynamic equations were found based only on the assumption that the endoreversible heat engine has a stationary state and from a Taylor expansion of the input/output heat fluxes in the first order. A relationship between relaxation times and total operation time was obtained. It was shown that the trajectories that lead the system back to the stationary state require much shorter times than the cycle time, allowing the system to work under continuous cyclic operation.
Mapping of the relaxation trajectories into the energetic space allowed for an analysis of the performance consequences of stability. The degree of symmetry of the conductance ratio in the input/output heat exchange is the main valuable factor for stability dynamics. For small σ hc values, the thermodynamic trajectories improve the efficiency and power values with a decrease in entropy generation, evolving near the endoreversible behavior. For larger σ hc values, the first part of the thermodynamic trajectories improve the efficiency and power values with a decrease in entropy generation, evolving toward the endoreversible behavior and the Pareto front. In between these two situations, i.e., for equal conductance values, (σ hc = 1) trajectories with decreasing efficiencies and increasing entropy generation can be found, with a preference for evolving near the Pareto front.
In summation, the stability of the irreversible Carnot-like heat engine exhibits two interesting behaviors in which the fluctuations around the stationary state (due to external perturbations) would likely maintain the system in an optimum state, or produce self-optimization induced by the stability. A biased control of the operation parameters could result in an economic saving, as energy is needed for the fine-tuning of parameter control. Finally, in this work the thermodynamic functions that were selected as the most relevant for this heat engine model were η, P, andṠ. It is very likely that by analyzing different stationary states, such as those predicted by the compromise Ecological [41] or Omega [42] operation performances, the role of the Pareto front would be more evident, as occurred in the case of the low-dissipation heat engine and refrigerator [1,2].
Another valuable remark is that the total entropy generation defines a time scale for both the operation and relaxation times. Finally, this analysis laid down the grounds to analyze heat devices dependent on working fluid properties, where the endoreversible hypothesis plays a relevant role. Especially when the thermalization mechanisms are fast enough compared to the cycle time in agreement with local equilibrium. This allows to incorporate, in a straightforward way, stochastic-type perturbations into the analysis through additive noise in Equation (14).