*Article* **Averaged Optimization and Finite-Time Thermodynamics**

## **Anatoly Tsirlin 1,\* and Ivan Sukin 2,\***

Ailamazyan Program Systems Institute of Russian Academy of Sciences, Petra Pervogo st., 4a, Veskovo, Yaroslavl oblast 152021, Russia

**\*** Correspondence: tsirlin@sarc.botik.ru (A.T.); ivsukin@gmail.com (I.S.); Tel.: +7-915-973-44-42 (I.S.)

Received: 1 July 2020; Accepted: 18 August 2020; Published: 20 August 2020

**Abstract:** The paper considers typical extremum problems that contain mean values of control variables or some functions of these variables. Relationships between such problems and cyclic modes of dynamical systems are explained and optimality conditions for these modes are found. The paper shows how these problems are linked to the field of finite-time thermodynamics.

**Keywords:** averaged; optimization; thermodynamics; heat transfer; cyclic mode

## **1. Introduction**

Averaging plays a very important role in optimization problems applied to engineering. Here are a few examples:

1. Assume that there is a bound on the source flow *gs* and this constraint lowers the possible optimal value of the production. If we introduce some buffer (container) in such a way that the source flow is its feed, we can raise the possible value of the feed flow to the actual process without violating this constraint. Actual values of the feed flow will oscillate between values that are greater than and less than *gs*. Only the mean value of the feed flow will be bounded in this case. Using this approach we can replace the strict constraint on *gs* by the averaged one. If we use such a buffer to store the product flow, we can maximize not this flow itself, but its mean value (Figure 1).

**Figure 1.** Flowsheet of a simple process with averaging of both source and product flows.

Let us assume that the relationship between production rate *g* and consumption *q* has the form presented at Figure 2.

**Figure 2.** Relationship between production and consumption. Effect of averaging.

With the help of buffers this relationship could be improved on interval from 0 to *q*1. The process must operate with consumption *q*<sup>1</sup> during some fraction of time and with zero consumption during the remaining time. Relationship between average production and average consumption is represented by dashed slope line at Figure 2. This is the way pumps of water towers operate.


Averaged problems arise in finite-time thermodynamics for two main reasons:


$$\frac{dX}{dt} = F(y) \tag{1}$$

The right hand side of (1) does not contain *X* and this means that the increase in extensive variables during some given amount of time depends only on the mean value of *F*. It does not depend on the order in which intensive variables have different values, if the mean value of *F* remains constant. Equations such as (1) are called Lyapunov-type. They allow us to formulate the problem of optimal control for thermodynamic systems in averaged form.

Below we will consider some of these problems. The last section contains applications of methods developed in the paper to finite-time thermodynamics.

We consider dynamical systems characterized by a finite number of variables.

By a *steady-state mode* of a system we mean a mode such that, for every variable *yν*(*t*) characterizing the system, there exists a period *T<sup>ν</sup>* such that the average value of *yν*(*t*) over this period is constant in time. Formally, this can be written as [1]

$$\frac{1}{T\_{\nu}} \int\_{t-T\_{\nu}}^{t} y\_{\nu}(\tau) d\tau = \overline{y}\_{\nu}. \tag{2}$$

Clearly, static modes, under which *yν*(*t*) are constant for all *ν*, satisfy this definition.

Another, more general, subclass of steady-state modes is formed by modes for which there exists a period *T* that is a multiple of all periods *Tν*. Such modes are called *cyclic* ones.

There are also steady modes for which there does not exist a common period *T* for all variables *yν*(*t*). This corresponds to the case where the ratio of at least two periods *T<sup>ν</sup>* and *T<sup>μ</sup>* is irrational. Such modes are called *quasi-cyclic* steady-state modes.

If the system is affected by external factors represented by stationary random processes and the mean values of the variables characterizing the system tend to some limits as the period *T* of averaging increases, then the steady mode is said to be *stochastic*.

A switch to a non-static steady-state mode may be caused either by the absence of a static mode admitted by the operating conditions of the system or by the fact that the efficiency of the system in a static mode is lower than in other modes.

Consider some examples.


In the rest of the paper, we mainly consider cyclic steady-state modes, among which two limit classes are distinguished. The first class includes modes in which each of the periods *Tν* significantly exceeds the time of relaxation processes in the system. Moreover, each of the static steady states is assumed to be stable. In this case, we can neglect the dynamics of the system and assume that under variation of the mode variables, the state variables change in accordance with the static characteristics. Such modes are said to be *quasi-static*.

The second class is formed by *sliding* steady-state modes, in which all or some of the control variables vary with frequency so high that, due to the inertia of the object, the state variables remain virtually constant, and their values depend only on the averaged influence of the control variables.

Although static modes are a special case of cyclic modes, below by cyclic modes we mean modes under which at least one variable of the process changes periodically in time. A cyclic mode is said to be efficient if the passage to this mode improves the efficiency of the process in comparison with the static mode.

3. Cyclic modes are typical of systems with no admissible static modes. Often a system has no static modes if the set *V* of admissible values of variables is non-convex; e.g., this set may include only discrete values. This is so, for example, in a heat engine in which a working fluid contacts a heat source whose temperature can take only two values, *T*<sup>+</sup> (a hot source) and *T*<sup>−</sup> (a cold source), and the average power over a cycle is required to be maximal under certain constraints.

Cyclic processes may be organized not only in time; variables may also depend on a spatial coordinate. In this case, the parameters of the system are constant in each section of the apparatus and vary periodically from section to section.

When passing from a static mode to a cyclic mode, one needs to replace the objective function by its mean value over the cycle and to replace all or some constraints imposed at each moment of time by averaged constraints. Thus, this passage involves an operation of averaging. Before passing to a cyclic mode, one must answer the following questions:


It is desirable to answer questions 1—3 without solving problem 4, which is rather difficult in most cases.

Usually the problem of choosing an optimal static mode of an apparatus reduces to the problem of finding the extremal value of the objective function under certain equality and inequality constraints on the variables, i.e., to a non-linear programming problem. The transition to cyclic modes extends the set of possible solutions and, depending on a particular setting, leads either to an averaged non-linear programming problem or to a variational control problem. If the problem has an optimal static mode, then we refer to this problem as the initial problem.

In the rest of the paper we consider various methods for constructing problems with larger sets of admissible solutions as compared with the initial problem; we show that there are relationships between such extended problems, which allows one to estimate solutions and values of some of them by solving others.

## **2. Averaged Optimization Problems and Their Optimality Conditions**

In this section, we consider various methods for introducing averaging into a non-linear programming (NLP) problem and obtain optimality conditions for averaged problems. To obtain these conditions, we use a trick based on reducing any averaged problem to a canonical form and deriving necessary optimality conditions for a particular problem from those for a general problem.

## *2.1. Averaging of Functions Included in the Formulation of an Optimization Problem*

Consider an initial NLP problem [2] in the form

$$f\_0(\mathbf{x}) \to \max \left/ f\_i(\mathbf{x}) = 0, \quad i = \overline{1, m}, \quad \mathbf{x} \in V\_\mathbf{x}. \tag{3}$$

On the set *Vx* , we define a probability measure *p*(*x*) such that

$$\int\_{V\_x} p(\mathbf{x})d\mathbf{x} = 1, \quad p(\mathbf{x}) \ge 0. \tag{4}$$

The average value of the function *f*(*x*) on the interval [0, *τ*] can be calculated as follows:

$$\overline{f(\mathbf{x})} = \frac{1}{\tau} \int\_0^\tau f(\mathbf{x}(t))dt = \int\_{V\_x} f(\mathbf{x})p(\mathbf{x})d\mathbf{x}.\tag{5}$$

Let us assume that *x* varies with time or one solves the problem (3) and maximizes the mean value of *f*0, but not the value of this function itself. If functions *fi* vanish on average, then we will arrive at a problem of the form

$$
\overline{f\_0(\mathbf{x})} \to \max \left/ \overline{f\_i(\mathbf{x})} = 0, \quad i = \overline{1, m}. \tag{6}
$$

A sought solution of problem (6) is a measure *p*(*x*) on *Vx* rather than a vector *x*. The variable *x* is called a *randomized* one, and *p*(*x*) is called a *generalized* solution. Following A.D. Ioffe and V.M. Tikhomirov [3], we call the value of the objective functional at the optimal solution as *the value of a problem*.

## *2.2. Convex Hulls—Carathéodory's Theorem*

The notion of *convexity* is very important for optimization problems.


Carathéodory's theorem is the most important theorem of convex analysis and geometry. It states that coordinates of every point of the convex hull of the set *<sup>V</sup>* <sup>⊂</sup> <sup>R</sup>*<sup>n</sup>* could be calculated as the weighted arithmetic mean of some points of *V* and the maximum necessary number of these points is no more than *n* + 1. The beautiful exposition of this theorem is given in [4].

## *2.3. Optimal Distribution in An Averaged NLP Problem*

Let us take some *x*<sup>0</sup> ∈ *Vx*. If *p*(*x*) = *δ*(*x* − *x*0), then problem (6) coincides with the initial problem. If the set of admissible solutions of a problem includes the set of admissible solution of the initial NLP problem and the optimality criteria in both problems coincide on the set of admissible solutions of the NLP problem, then the former problem is called an extension of the NLP problem.

First, consider the special form of problem (6) with *fi*(*x*) = *xi*:

$$
\overline{f\_0(\mathbf{x})} \to \max \left/ \overline{\mathbf{x}}\_i = 0, \quad i = \overline{1, n}. \tag{7}
$$

The value of the problem (7) is equal to the ordinate of the convex hull of the function *f*0(*x*) on the set *Vx* at the point *x* = 0. According to Carathéodory's theorem, constructing any ordinate of the convex hull of a function of n variables requires averaging at most *n* + 1 ordinates of the function *f*0(*x*); therefore, we can rewrite problem (7) in the form

$$\sum\_{\nu=0}^{n} \gamma\_{\nu} f\_0(\mathbf{x}^{\nu}) \to \max \left/ \sum\_{\nu=0} \gamma\_{\nu} \mathbf{x}\_i^{\nu} = 0, \quad i = \overline{1, n}, \quad \gamma\_{\nu} \ge 0, \quad \sum\_{\nu=0}^{n} \gamma\_{\nu} = 1. \tag{8}$$

Let us return to problem (6) and try to reduce it to simple calculation. We need to calculate the ordinate of a convex hull of the given function. Please note that problem (6) can be solved in two stages. At the first stage, we find the maximum of the function *f*0(*x*) subject to the constraint *f*(*x*) = *C*, where *C* takes all values for which the level surface *f*(*x*) = *C* intersects *Vx*. The problem

$$f\_0(\mathbf{x}) \to \max \left/ f\_i(\mathbf{x}) = \mathbb{C}, \quad i = \overline{1, m}, \quad \mathbf{x} \in V\_x \tag{9}$$

is a non-linear programming problem. Solving (9), we obtain a set of conditionally optimal solutions *x*∗(*C*) and the corresponding values of the *reachability function f* ∗ <sup>0</sup> (*C*) = *f*0(*x*∗(*C*)) of the non-linear programming problem.

The following assertion holds: *The optimal distribution p*∗(*x*) *in problem (6) is concentrated at the points x*∗(*C*). In other words, *one needs to average only over conditionally optimal values of f*0.

## *2.4. Necessary Conditions of Optimality—Kuhn-Tucker Theorem*

The Kuhn-Tucker theorem generalizes Lagrange multipliers method to problems with inequality constraints:

$$f\_0(\mathbf{x}) \to \max\_{\mathbf{x}} \left/ f\_i(\mathbf{x}) = 0, \quad \mathbf{q}\_j(\mathbf{x}) \ge 0, \quad i = \overline{1, k}, \quad j = \overline{k+1, m} \right. \tag{10}$$

where all functions are smooth.

The theorem states that *there is nonzero vector of Lagrange multipliers with components λi, μ<sup>j</sup>* ≤ 0 *such that Lagrange function*

$$L = \lambda\_0 f\_0(\mathbf{x}) + \sum\_i \lambda\_i f\_i(\mathbf{x}) + \sum\_j \mu\_j \varphi\_j(\mathbf{x}) = R(\lambda, \mathbf{x}) + \sum\_j \mu\_j \varphi\_j(\mathbf{x}) \tag{11}$$

*is stationary on the optimal solution of the problem (10). The multiplier λ*<sup>0</sup> *could equal to zero or one. In the former case the solution is called degenerate.*

It follows from this theorem that when *ϕj*(*x*) = *xj* we have inequality *<sup>∂</sup><sup>R</sup> <sup>∂</sup>xj* ≤ 0 for the optimal solution. More detailed explanation could be found in [5].

## *2.5. Reduction to an Ordinary NLP Problem*

The above considerations allow us to formulate the second stage in solving problem (6). This stage is the maximization of the average value of the function *f* ∗ <sup>0</sup> (*C*) with the constraint that the vector *C* has zero mean, i.e.,

$$\sqrt[4]{f\_0^\*(\mathbb{C})} \to \max \sqrt{\mathbb{C}\_i} = 0, \quad i = \overline{1, m}, \quad \mathbb{C}\_i \in V\_{\mathbb{C}}.\tag{12}$$

This problem is similar to the problem (7). Its value, and hence the value of problem (6), is equal to the ordinate of the convex hull of the reachability function *f* ∗ <sup>0</sup> (*C*) at *C* = 0:

$$\sup\_{\mathbf{x}\in D} \overline{f\_0(\mathbf{x})} = \sup \overline{f\_0^\*(\mathbf{C})} \; \bigg/ \sqrt{\mathbf{C}} = \mathbf{0}, \quad \mathbf{C} \in V\_{\mathbf{C}}.\tag{13}$$

Since the vector *C* is *m*-dimensional, the number of base points *C<sup>ν</sup>* in problem (9) is at most *m* + 1. Thus, the distribution *p*(*C*) in problem (9) can be sought in the form

$$p(\mathbb{C}) = \sum\_{\nu=0}^{m} \gamma\_{\nu} \delta(\mathbb{C} - \mathbb{C}^{m}). \tag{14}$$

Since each of the base values *C<sup>ν</sup>* corresponds to a conditionally optimal solution *x*∗(*Cν*), the optimal distribution *p*(*x*) is also concentrated at no more than *m* + 1 points:

$$p(\mathbf{x}) = \sum\_{\nu=0}^{m} \gamma\_{\nu} \delta(\mathbf{x} - \mathbf{x}^{\nu}). \tag{15}$$

Substituting the distribution (15) into the expressions for *f*0(*x*) and *fi*(*x*), we reduce problem (6) to the form

$$\begin{aligned} I &= \sum\_{\nu=0}^{m} \gamma\_{\nu} f\_0(\mathbf{x}^{\nu}) \to \max \left/ \sum\_{\nu=0}^{m} \gamma\_{\nu} f\_i(\mathbf{x}^{\nu}) = 0, \\\ \mathbf{i} &= \overline{1, m}, \quad \mathbf{x}^{\nu} \in V\_{\mathbf{x}\prime} \quad \gamma\_{\nu} \ge 0, \quad \sum\_{\nu=0}^{m} \gamma\_{\nu m} = 1. \end{aligned} \tag{16}$$

Thus, we have reduced the problem to an ordinary NLP problem whose variables are the base values *x<sup>ν</sup>* of the vector *x* and the weight factors *γν*.

*2.6. Relationship between Averaged NLP Problem and the Lagrangian Function of the NLP Problem without Averaging*

The Lagrangian function

$$\begin{array}{lcl}\overline{R} &= \sum\_{\nu=0}^{m} \gamma\_{\nu} f\_{0}(\mathbf{x}^{\nu}) + \sum\_{i=1}^{m} \lambda\_{i} \sum\_{\nu=0}^{m} \gamma\_{\nu} f\_{i}(\mathbf{x}^{\nu}) + \Lambda \left(1 - \sum\_{\nu=0}^{m} \gamma\_{\nu}\right) = \\ & \sum\_{\nu=0}^{m} \gamma\_{\nu} \left[f\_{0}(\mathbf{x}^{\nu}) + \sum\_{i=1}^{m} \lambda\_{i} f\_{i}(\mathbf{x}^{\nu}) - \Lambda\right] + \Lambda \end{array} \tag{17}$$

of problem (16) is related to the Lagrangian function

$$R = f\_0(\mathbf{x}) + \sum\_{i=1}^{m} \lambda\_i f\_i(\mathbf{x}) \tag{18}$$

of the initial NLP problem by

$$\overline{R} = \sum\_{\nu=0}^{m} \gamma\_{\nu} (R(\mathbf{x}^{\nu}, \lambda) - \Lambda) + \Lambda. \tag{19}$$

Since ∑ *ν γν* = 1, the Lagrangian function of the averaged problem equals the average value of the Lagrangian function of the initial problem over all base values *xν*. Some of the weight factors *γν* may vanish; then the number of base points is less than *m* + 1.

Let us find conditions that must hold for those *x<sup>ν</sup>* that have non-zero weights in (19). For this purpose, we apply the Kuhn–Tucker theorem and write the optimality conditions for problem (16) with respect to the variables *γν*:

$$
\frac{\partial \overline{\mathcal{R}}}{\partial \gamma\_{\nu}} \delta \gamma\_{\nu} \le 0. \tag{20}
$$

Since *γν* are bounded only from below (*γν* ≥ 0), it follows that *δγν* ≥ 0; therefore,

$$\frac{\partial \overline{R}}{\partial \gamma\_{\nu}} = R(\mathbf{x}^{\nu}) - \Lambda \le 0,\tag{21}$$

or *<sup>R</sup>*(*xν*) <sup>≥</sup> <sup>Λ</sup>. If *<sup>γ</sup>*<sup>∗</sup> *<sup>ν</sup>* > 0, then *δγν* may be of any sign, and so inequality (21) transforms into the equality

$$R(\mathbf{x}^{\nu}) = \Lambda. \tag{22}$$

Thus, *for all x<sup>ν</sup> involved in the averaged problem with non-zero weights, the Lagrangian function R of the initial non-linear programming problem attains an absolute maximum*. Of course, this maximum is the same for all *xν*.

The requirements that the function *R* must take the same value at all points *xν*<sup>∗</sup> and that this value must be maximum give equations for the variables to be found. Thus, applying Kuhn–Tucker theorem the problem (6), we obtain the vector of Lagrange multipliers *λ* for which the function *R* attains an absolute maximum with respect to the variables *<sup>x</sup><sup>ν</sup>* <sup>∈</sup> *Vx* and *γν* <sup>∈</sup> *<sup>V</sup><sup>γ</sup>* at an element of the set *<sup>D</sup>* of admissible solutions to problem (6), and these multipliers *λ* satisfy the condition

$$\overline{\mathcal{R}}(\lambda^\*, \gamma\_{\nu^\*}^\* \mathbf{x}^{\nu \*}) = \inf\_{\lambda \in V\_{\lambda}} \sup\_{\gamma\_{\nu^\*} \mathbf{x}^{\nu}} \overline{\mathcal{R}}(\lambda, \gamma\_{\nu^\*} \mathbf{x}^{\nu}) = \inf\_{\lambda \in V\_{\lambda}} \sup\_{\mathbf{x} \in V\_{\lambda}} \mathcal{R}(\lambda, \mathbf{x}). \tag{23}$$

Thus, when the attainability function *f* ∗ <sup>0</sup> (*C*) coincides with its convex hull at *C* = 0, the transition to the averaged problem is not efficient (the values of the NLP problem and problem (6) coincide). By virtue of (23), we can look for the value of the averaged problem in the form inf *λ*∈*V<sup>λ</sup>* sup *x*∈*Vx R*(*λ*, *x*). If the

extended problem is inefficient, then we say that it is equivalent to the initial problem.

In the general case, the dimension of the vector of unknown variables and the computational complexity of problem (6) are much greater than those for the NLP problem. However, in many cases, we are interested not in the solution but in the value of the averaged problem, which shows the gain obtained by transition to the averaged setting. Some methods for estimating the value of problem (6) from above and below were proposed in [6].

#### *2.7. Other Forms of Averaged Extensions of the NLP Problem*

Problem (6) is not the only possible extension of the NLP problem by averaging. The optimality criteria, relations, and constraints in real-life problems often include the mean values of variables *x* rather than the variables themselves. For example, the performance of a distillation column

is characterized by the mean not current composition of output flows, because these flows are accumulated in some containers or apparatuses at the exit of the column (or attached to the column). Below, we describe several possible modifications of the averaged extension [7].

1. *Problem of maximizing a function of the mean value of the argument.* When *D* is the set of admissible solutions of the initial NLP problem, i.e., *D* is defined by the condition *f*(*x*) = 0, and *x* is the mean value of the vector *x* on the set *D*, we have:

$$f\_0(\mathbf{x}) \to \sup \slash p(\mathbf{x}) = 0 \quad \forall \mathbf{x} \notin D. \tag{24}$$

Since the set of values *x* satisfying this condition is the convex hull of *D*, problem (24) is equivalent to the NLP problem on the convex hull of *D*:

$$f\_0(\mathbf{x}) \to \sup \mid \mathbf{x} \in \mathbf{Co} \, D. \tag{25}$$

2. *Problem of maximizing the mean value of a function under constraints imposed on the mean value of the argument:*

$$\begin{array}{c}\overline{f\_0(\mathbf{x})} \rightarrow \sup \Big/ f\_i(\overline{\mathbf{x}}) = 0, \quad i = \overline{1, m} \end{array} \tag{26}$$

or, in more detail,

$$\int\_{V\_x} f\_0(\mathbf{x}) p(\mathbf{x}) d\mathbf{x} \to \sup\_{p(\mathbf{x})} \left\langle f\_i \left( \int\_{V\_x} \mathbf{x} \cdot p(\mathbf{x}) d\mathbf{x} \right) = 0, \quad i = \overline{1, m}. \tag{27}$$

3. *Problem of maximizing a function of the mean value of x under averaged constraints:*

$$f\_0(\overline{\mathfrak{x}}) \to \sup \left/ \overline{f\_i(\mathfrak{x})} = 0, \quad i = \overline{1, m}. \tag{28}$$

Each of the above problems is an extension of the non-linear programming problem, and the solutions of these problems are distributions *p*(*x*).

*Averaged problems with two types of variables*. An NLP problem can be extended only with respect to some components of the solution rather than with respect to the whole solution. In practice, this situation occurs when the problem is solved repeatedly and some components (we denote them by *x*) can vary from one solution to another, while the remaining components must be chosen only once and then fixed. We denote the latter group of variables by *y*. For example, *x* may be the operating conditions of the process (such as flow, pressure, temperature, etc.) and *y* may be the design parameters of an apparatus.

If we denote

$$\overline{f(y,\mathbf{x})}^{\mathbf{x}} = \int\_{\mathbb{C}\_x} f(y,\mathbf{x}) p(\mathbf{x}) d\mathbf{x},\tag{29}$$

a problem in which averaging is performed over only part of variables has the form

$$\overline{f\_0(y,\mathbf{x})}^x \to \sup \left/ \overline{f\_i(y,\mathbf{x})}^x = 0, \quad i = \overline{1,m}.\tag{30}$$

One need to find the vector *y* and distribution *p*(*x*) in (30).

For each fixed *y*, this problem coincides with the usual setting of problem (6). If we separate the randomized variables *<sup>x</sup>* <sup>∈</sup> *<sup>E</sup><sup>r</sup>* and the deterministic variables *<sup>y</sup>* <sup>∈</sup> *<sup>E</sup><sup>s</sup>* in the Lagrangian function *<sup>R</sup>* of the initial NLP problem, then we can write optimality conditions with respect to *x* by analogy with problem (6) in the form (see (23))

$$R(\lambda, \gamma\_{\nu^{\prime}}^{\*} y, \mathbf{x}^{\nu \*}) = \sup\_{\mathbf{x} \in V\_{\mathbf{x}}} R(\lambda, y, \mathbf{x}), \quad \nu = \overline{0, m}. \tag{31}$$

In this case, if we denote the admissible set of (30) as *Dx*(*y*) *x* , for each *y* ∈ *Vy*, there exist *λ*(*y*) such that

$$\inf\_{\lambda} \sup\_{x \in V\_x} R(\lambda, y, x) = \sup\_{x \in \overline{D\_x(y)}^x} \overline{f\_0(y, x)}.\tag{32}$$

The Lagrangian function attains an absolute maximum at the base values of *x*.

At the same time, for a fixed function *p*(*x*), problem (30) becomes a usual non-linear programming problem with respect to the variables *y*. The Kuhn–Tucker conditions hold for this problem, which include in this case the complementary slackness conditions as well as the requirement that the function *R*(*λ*, *γν*, *y*, *xν*) be stationary with respect to *y*, which in turn, leads to the equations

$$\frac{\partial}{\partial y\_j} \left[ \sum\_{\nu=0}^m \gamma\_\nu \mathcal{R}(\lambda, y, \mathbf{x}^\nu) \right] = 0, \quad j = \overline{1, s}. \tag{33}$$

where *R* is the Lagrangian function for the NP problem.

Averaged problems with two types of variables are in a sense close to optimal control problems, and optimality conditions for such problems are close to the Pontryagin maximum principle.

## *2.8. The Algorithm for Obtaining Optimality Conditions in Averaged Problems*

By *an averaged problem of static optimization* we mean any NLP problem in which either functions or variables are averaged with respect to all or part of the variables.

As shown above, the settings of averaged problems are very diverse. The reason for this is that a problem may contain both the mean values of functions and functions of the mean values of variables. Moreover, averaging may involve only part of the variables. Under these conditions, it is inexpedient to study each possible setting of an averaged problem. It is significantly more convenient to obtain optimality conditions for some canonical form of an averaged problem and apply them to each particular problem after having reduced the latter to this canonical form [8].

Before obtaining optimality conditions, we must answer the following two questions:


The necessary optimality conditions given below yield an affirmative answer to the first question and allow one to determine the limit number of base points.

Let *y* denote the vector of deterministic variables, and let *x* be the vector of randomized variables. For the former, we must find an optimal value, and for the latter, an optimal measure. The canonical form of the averaged problem is

$$F\_0(\overline{f(\mathbf{x}, y)}, y, \mathbf{x}) \to \max \tag{34}$$

under the constraints

$$\begin{array}{ll} F\_{\boldsymbol{\nu}}(\overline{f(\boldsymbol{x},\boldsymbol{y})},\overline{\boldsymbol{\rho}(\boldsymbol{x},\boldsymbol{y})},\overline{\boldsymbol{\pi}})=0, & \boldsymbol{\nu}=\overline{1,r},\\ F\_{\boldsymbol{\nu}}(\overline{f(\boldsymbol{x},\boldsymbol{y})},\overline{\boldsymbol{\rho}(\boldsymbol{x},\boldsymbol{y})},\overline{\boldsymbol{\pi}})\geq 0, & \boldsymbol{\nu}=\overline{r+1,m}.\end{array} \tag{35}$$

Here the bar over the symbol of a function denotes averaging over the set *Vx* of randomized variables *x*, which is assumed to be compact.

Suppose that the vector *x* has dimension *k* and the vector function *f* has dimension *n*. The function *F* is assumed to be continuously differentiable with respect to all its variables, and *f* and *ϕ* are continuous in *x* and continuously differentiable in *y*.

In [8], one of the authors (*A.T.*) proved that *the optimal measure p*∗(*x*) *on the set of randomized variables is concentrated at no more than L* + 1 *base points, where L* = *n* + *k*. Thus,

$$p^\*(\mathbf{x}) = \sum\_{l=0}^L \gamma\_l \delta(\mathbf{x} - \mathbf{x}^l), \quad \gamma\_l \ge 0, \quad \sum\_{l=0}^L \gamma\_l = 1. \tag{36}$$

Therefore, for the optimal solution, we have

$$\overline{f^\*(\mathbf{x}, \mathbf{y})} = \sum\_{l=0}^L \gamma\_l f(\mathbf{x}^l, \mathbf{y}), \quad \overline{\mathbf{x}} = \sum\_{l=0}^L \gamma\_l \mathbf{x}^l,\tag{37}$$

and constraints (35) take the form

$$\begin{array}{l} F\_{\nu}(\overline{f}, \overline{\varrho(\mathbf{x}^{l}, y)}, \overline{\mathbf{x}}) = 0, \quad \nu = \overline{1, r}\_{\prime}, \\ F\_{\nu}(\overline{f}, \overline{\varrho(\mathbf{x}^{l}, y)}, \overline{\mathbf{x}}) \ge 0, \quad \nu = \overline{r + 1}\_{\prime}, \end{array} \tag{38}$$

for all values of *x<sup>l</sup>* .

These expressions turn problem (34), (35) into an ordinary NLP problem with respect to *γl*, *y* and *xl* . The Kuhn–Tucker conditions reduce to the following: the Lagrangian function

$$R = F\_0(\overline{f}, y, \overline{\pi}) + \sum\_{\nu=1}^{m} \lambda\_\nu F\_\nu(\overline{f}, \varphi(x^I, y), \overline{\pi}) \tag{39}$$

of this problem is stationary with respect to *x<sup>l</sup>* and *y* and is unimprovable with respect to *γ<sup>l</sup>* (we assume the solution is non-degenerate, so *λ*<sup>0</sup> = 1). To write down the optimality conditions, we introduce the notation

$$a\_{\!\!\!\!\!} = \frac{\partial \mathcal{R}}{\partial \overline{f}\_{\!\!\!\!}}, \quad \beta\_{\!\!\!\!\/} = \frac{\partial \mathcal{R}}{\partial x\_{\!\!\!\!\/}}, \quad r\_{\mu \!\!\!\/]} = \frac{\partial \mathcal{R}}{\partial q\_{\mu}}(x^{\!\!\!\/]} \mathbf{x}^{\!\!\!\!/} \mathbf{y}). \tag{40}$$

Using this notation, we can write the condition that *R* is unimprovable with respect to *γ<sup>l</sup>* as follows: the expression

$$\mathcal{C}(\mathbf{x}) = \sum\_{j} a\_{j} f\_{j}(\mathbf{x}, \mathbf{y}^{\*}) + \sum\_{i} \beta\_{i} \mathbf{x}\_{i} \tag{41}$$

attains its maximum with respect to *<sup>x</sup>* <sup>∈</sup> *Vx* at the points *<sup>x</sup><sup>l</sup>* , so that

$$\mathbf{x}^{l\*} = \underset{\mathbf{x}}{\arg\max} \mathcal{C}(\mathbf{x}), \quad l = \overline{1, L}; \tag{42}$$

the condition that *R* is stationary with respect to *y* has the form

$$\nabla\_y \left[ \sum\_j a\_j \overline{f\_j(\mathbf{x}, y)} + F\_0(\overline{f}, \mathbf{x}, y) + \sum\_{\mu, l} r\_{\mu l} q\_{\mu}(\mathbf{x}^l, y) \right] = 0. \tag{43}$$

The maximality of *C*(*x*), together with equations (42), constraints (35), and the complementary slackness conditions

$$\sum\_{\nu=r+1}^{m} \lambda\_{\nu} F\_{\nu}(\overline{f}^\*, \overline{\rho}^\*, \mathbf{x}^\*) = 0, \quad \lambda\_{\nu} \ge 0, \quad \nu = \overline{r+1, m} \tag{44}$$

allows one to find a solution *γ*∗ *<sup>l</sup>* , *<sup>y</sup>*∗, *<sup>x</sup><sup>l</sup>* .

When formulating a specific averaged problem, one


For example, in problem (26), we have

$$F\_0 = \overline{f\_0(\mathbf{x})}, \quad F\_\nu = f\_\nu(\overline{\mathbf{x}}), \quad \nu = \overline{1, m}. \tag{45}$$

The number *L* equals *k*, and

$$R = \lambda\_0 \overline{f\_0(\mathbf{x})} + \sum\_{\nu=1}^{m} \lambda\_\nu f\_\nu(\overline{\mathbf{x}}). \tag{46}$$

In (42), we have *a*<sup>0</sup> = *λ*<sup>0</sup> = 1, *a<sup>ν</sup>* = 0 for *ν* > 0 and

$$\beta\_i = \sum\_{\nu=1}^m \lambda\_\nu \left( \frac{\partial f(\overline{\mathbf{x}})}{\partial \overline{\mathbf{x}\_i}} \right)\_{\mathbf{x}'} \quad i = \overline{1,k}. \tag{47}$$

At the base points *x<sup>l</sup>* , the number of which does not exceed *k* + 1, the expression

$$\mathcal{C}(\mathbf{x}) = f\_0(\mathbf{x}) + \sum\_{i=1}^k \mathbf{x}\_i \sum\_{\nu=1}^m \lambda\_\nu \left( \frac{\partial f(\overline{\mathbf{x}})}{\partial \overline{\mathbf{x}}\_i} \right)\_{\mathbf{x}^\*} \tag{48}$$

attains its maximum, and conditions (35) hold, which have the form

$$f\_{\nu}\left(\sum\_{l=0}^{k}\gamma\_{l}\mathbf{x}^{l}\right)=0,\quad\nu=\overline{1,m}.\tag{49}$$

#### **3. Non-Stationary Problems of Averaged Optimization**

Consider an extremal problem of the form

$$\overline{f}\_0 = \frac{1}{\pi} \int\_0^\pi f \mathbf{o}(J(t), u(t)) dt \to \max\_{u} \tag{50}$$

subject to the constraints

$$\overline{f}\_{\nu} = \frac{1}{\pi} \int\_{0}^{\overline{\tau}} f\_{\nu}(f(t), u(t)) dt = 0, \quad \nu = \overline{1, n}, \tag{51}$$

where the functions *<sup>f</sup><sup>ν</sup>* : <sup>R</sup>*k*<sup>1</sup> <sup>×</sup> <sup>R</sup>*k*<sup>2</sup> <sup>→</sup> <sup>R</sup>, *<sup>ν</sup>* <sup>=</sup> 0, *<sup>n</sup>*, are continuous in *<sup>J</sup>* and *<sup>u</sup>*, *<sup>u</sup>* <sup>∈</sup> *Vu* <sup>⊂</sup> <sup>R</sup>*k*<sup>1</sup> is a measurable function, the set *Vu* is compact, and *<sup>J</sup>*(*t*) <sup>∈</sup> *VJ* ⊂ R*k*<sup>2</sup> is a given measurable function of time. With *J*(*t*) we can associate a probability measure (distribution) *p*(*J*). If *J*(*t*) takes a value *J<sup>k</sup>* on a part of the interval (0, *<sup>τ</sup>*) of relative length *<sup>α</sup>k*, then *<sup>p</sup>*(*J*) contains a term of the form *<sup>α</sup>kδ*(*<sup>J</sup>* <sup>−</sup> *<sup>J</sup>k*). The length of the interval (0, *τ*) may tend to infinity, and *J*(*t*) may be a stationary random process with distribution *p*(*J*).

The distribution *p*(*J*) can be written in the form

$$p(f) = \overline{p}(f) + \sum\_{k} a\_k \delta(f - f^k). \tag{52}$$

For problem (50), (51), let *ατ* be the length of the part of (0, *τ*) on which *J*(*t*) takes one of the constant values *Jk*; we have *α* = ∑ *k αk*. We refer to *ατ* as the total constancy interval of *J*(*t*). The remaining part (1 − *α*)*τ* is called the interval of variation of the parameter *J*.

**Theorem 1.** *Let u*∗(*t*) *be an optimal solution; then there exists a non-zero vector λ* = {*λ*0, ... , *λn*} *with λ*<sup>0</sup> ∈ {0, 1} *such that*

• *on the interval of variation of the parameter J*(*t*)

$$\mu^\*(f,\lambda) = \underset{u \in V\_H}{\text{arg}\max} \sum\_{\nu=0}^n \lambda\_\nu f\_\nu(f, u);\tag{53}$$

• *on the total constancy interval of J*(*t*)*, the optimal solution switches between at most n* + 1 *base values u<sup>j</sup> , and each of these values satisfies the condition*

$$\mu^j = \underset{u \in V\_u}{\text{arg}\max} \sum\_k a\_k \sum\_{\nu=0}^n \lambda\_\nu f\_\nu(f^k, \mu), \quad j = \overline{0, n};\tag{54}$$

• *the portions γ<sup>j</sup> of the constancy interval ατ on which u*∗(*t*) *takes the respective values u<sup>j</sup> satisfy the conditions*

$$\begin{aligned} \int\_{V} \overline{p}(I) f\_{\boldsymbol{\nu}}(I, \boldsymbol{u}^{\star}(I)) dI + \sum\_{j=0}^{n} \gamma\_{j} \sum\_{k} a\_{k} f\_{\boldsymbol{\nu}}(I^{k}, \boldsymbol{u}^{j}) &= 0, \quad \boldsymbol{\nu} = \overline{1, n}, \\ \sum\_{j=0}^{n} \gamma\_{j} &= 1, \quad \gamma\_{j} \ge 0; \end{aligned} \tag{55}$$

• *the vector of multipliers λν, ν* = 1, *n, is determined by the conditions*

$$\lambda^\* = \arg\min\_{\lambda} \left[ \int\_{\mathcal{V}} \overline{p}(I) \sum\_{\nu=0}^n \lambda\_{\nu} f\_{\nu}(f\_{\nu} \mu^\*(I\_{\nu} \lambda)) dI + \sum\_{j=0}^n \gamma\_j \sum\_{\nu=0}^n \lambda\_{\nu} \sum\_k a\_k f\_{\nu}(f^k, \mu^j(\lambda)) \right]. \tag{56}$$

Thus, on the constancy intervals, the optimal solution of a problem with non-stationary parameters coincides with the solution of an averaged mathematical programming problem, and on the interval of variation of the parameter, it varies as the solution of a problem with integral constraints. This theorem was proved in [9].

**Example 1.** *Consider the problem of maximizing the average power p of a heat engine in which the working fluid contacts a source of variable temperature T*0(*t*)*. This problem has the form*

$$\overline{p} = \frac{1}{\pi} \int\_0^\overline{q} (T\_0(t), T(t)) dt \to \max\_T \tag{57}$$

*subject to the constraint*

$$
\overline{\sigma} = \frac{1}{\tau} \int\_0^\tau \frac{q(T\_0(t), T(t))}{T(t)} dt = 0. \tag{58}
$$

*Here T*(*t*) *is the temperature of the working substance, q is the heat flux from the source to the working fluid, and σ is the mean rate of variation of the entropy of the working substance. A substantiation of the setting (56), (57) can be found in [10–12]. The optimality conditions (53) imply the following relation for the interval of variation of T*0(*t*)*:*

$$\frac{1}{T^2} \frac{q(T\_{0\prime}T)}{\partial q(T\_{0\prime}T)/\partial T} - \frac{1}{T} = \text{const.}\tag{59}$$

In particular, for the Newtonian law *q*(*T*0, *T*) = *β*(*T*<sup>0</sup> − *T*) of heating, (59) implies

$$T^\*(T\_0) = m\sqrt{T\_0} \tag{60}$$

where *m* is the constant equal to the mean value of the square root of the source temperature.

For example, suppose that *T*0(*t*) has a uniform distribution (for a regular function *T*0(*t*), this means that the source temperature depends linearly on time) and *T*<sup>02</sup> and *T*<sup>01</sup> are the maximal and minimal source temperatures, respectively. Then

$$T^\*(T\_0) = \frac{2(T\_{02}^{3/2} - T\_{01}^{3/2})}{3(T\_{02} - T\_{01})} \sqrt{T\_0}.\tag{61}$$

The maximum power is given by

$$\overline{p}\_{\text{max}} = \beta \left[ \frac{T\_{02} + T\_{01}}{2} - \frac{4}{9} \frac{(T\_{02}^{3/2} - T\_{01}^{3/2})}{T\_{02} - T\_{01}} \right]. \tag{62}$$

Thus, a heat engine with one source may have non-zero power if the variance of the source temperature is positive.

For some laws *q*(*T*0, *T*), the optimal temperature *T*∗(*t*) may switch between two base values on intervals of constancy of the parameter *T*0.

## **4. Estimation of the Performance of Cyclic Modes**

Suppose that the dynamics of a system is characterized by the differential equations

$$
\dot{\mathbf{x}}\_{\nu} = f\_{\nu}(\mathbf{x}, \boldsymbol{\mu}, \boldsymbol{a}), \quad \nu = \overline{1, m}\_{\prime} \tag{63}
$$

whose right-hand sides do not explicitly depend on *t*. Here, as in the preceding sections, *x* denotes the state variables, *u* are the control ones, and *a* denotes parameters to be optimized. As a rule, boundary conditions are not fixed for equations (63), but the state variables are required to vary periodically:

$$\mathbf{x}\_{\nu}(\tau) = \mathbf{x}\_{\nu}(0) \Rightarrow \int\_{0}^{\tau} f\_{\nu}(\mathbf{x}, u, a) dt = 0, \quad \nu = \overline{1, m}. \tag{64}$$

The performance averaged over the cycle plays the role of the optimality criterion for such a cyclic process and can be written in the form

$$I = \frac{1}{\pi} \int\_0^\pi f\_0(\mathbf{x}, u, a) dt \to \max. \tag{65}$$

The duration *τ* of each cycle is one of the components of the vector *a*; in the general case, it is not fixed. The parameters and controls are subject to constraints *a* ∈ *Va* and *u* ∈ *Vu*; in addition to the integral constraints (64), which follow from the periodicity of the process, the problem usually contains integral constraints determined by given mean rates of consumption of some resources (resource constraints):

$$J\_j = \int\_0^\tau \varphi\_j(\mathbf{x}, \mu, a)dt = 0, \quad j \in \overline{1, r}. \tag{66}$$

It is assumed that each of the functions determining the problem is continuous in all its variables and is continuously differentiable with respect to *x* and *a*.

**Optimality conditions**. Optimality conditions for problem (63)–(66) can be obtained by using the maximum principle [6]. Namely, if an optimal solution *x*∗, *a*∗, *u*∗ exists and is non-degenerate, then there exist a non-zero vector *λ* and a differentiable vector function *ψ*(*t*) such that the function

$$\mathcal{R} = \frac{1}{\pi} f\_0 + \sum\_{\nu} \left[ \dot{\psi}\_{\nu} \mathbf{x}\_{\nu} + (\psi\_{\nu} + \lambda\_{\nu}) f\_{\nu} \right] + \sum\_{j} \lambda\_{j} \boldsymbol{\uprho}\_{j}.\tag{67}$$

is stationary with respect to *x* and attains a maximum with respect to *u*, and the integral *S* of this function is locally unimprovable with respect to *a*. Thus,

$$\frac{\partial R}{\partial \mathbf{x}\_i} = 0 \quad \Rightarrow \quad \dot{\boldsymbol{\psi}}\_i = -\frac{\partial}{\partial \mathbf{x}\_i} \left\{ \frac{1}{\tau} f\_0 + \sum\_{\nu} (\psi\_{\nu} + \lambda\_{\nu}) f\_{\nu} + \sum\_{j} \lambda\_j \boldsymbol{\varphi}\_j \right\}. \tag{68}$$

Since the values *xν*(*τ*) and *xν*(0) are not fixed, it follows that *ψν*(*τ*) and *ψν*(0) vanish. Introducing the notation *ψ*˜*<sup>ν</sup>* = *ψν* + *λν* and taking into account the equality ˙ *ψ*˜*<sup>ν</sup>* = *ψ*˙ *<sup>ν</sup>*, we can rewrite condition (68) in the form

$$\hat{\boldsymbol{\psi}}\_{i} = -\frac{\partial}{\partial \mathbf{x}\_{i}} \left\{ \frac{1}{\pi} f\_{0} + \sum\_{\mathcal{V}} \psi\_{\mathcal{V}} f\_{\mathcal{V}} + \sum\_{j} \lambda\_{j} \boldsymbol{\varphi}\_{j} \right\} = -\frac{\partial}{\partial \mathbf{x}\_{i}} H. \tag{69}$$

For these equations, since *ψ*(0) and *ψ*(*τ*) vanish, the costate variables satisfy the periodicity conditions

$$
\tilde{\psi}\_{\nu}(0) = \tilde{\psi}\_{\nu}(\tau) \quad \Rightarrow \quad \int\_{0}^{\tau} \frac{\partial H}{\partial \mathbf{x}\_{\nu}} dt = 0, \quad \nu = \overline{1, m}. \tag{70}
$$

The conditions of maximality of *R* with respect to *u* have the form

$$u^\*(t) = \underset{u \in V\_{\mu}}{\text{arg}\max} \left\{ \frac{f\_0}{\tau} + \sum\_{\nu} \tilde{\psi}\_{\nu} f\_{\nu} + \sum\_{j} \lambda\_j \varphi\_j \right\}. \tag{71}$$

Finally, the optimality conditions with respect to each component *ak* of the vector *a*, including the duration *τ* of the cycle, yield the inequalities

$$\frac{\partial S}{\partial a\_k} \delta a\_k \le 0, \quad k = 1, 2, \dots \tag{72}$$

Here *δa* is the cone of variations of the vector *a* that are admissible with respect to the inclusion *a* ∈ *Va*.

Please note that the phase trajectory corresponding to an optimal cyclic process has no self-intersections [13].

#### **5. Estimation of the Efficiency of Transition to a Cyclic Process**

## *5.1. Conditions of Equivalence and Efficiency of a Cyclic Extension*

The optimal cyclic mode problem (63)–(66) (we refer to it as Problem C) is an extension of a non-linear programming problem. Indeed, imposing the additional constraints *x* = const and *u* = const on the solution of this problem, we obtain the following optimal static mode problem (Problem S):

$$d\_{\mathcal{S}} = f\_0(\mathbf{x}, \boldsymbol{\mu}, a) \to \max \left/ \begin{array}{c} f\_{\boldsymbol{\nu}}(\mathbf{x}, \boldsymbol{\mu}, a) = 0, \quad \boldsymbol{\varphi}\_{\boldsymbol{\bar{\jmath}}}(\mathbf{x}, \boldsymbol{\mu}, a) = 0 \\ \boldsymbol{\mu} \in V\_{\mathbf{u}\boldsymbol{\nu}} \quad a \in V\_{\mathbf{u}\boldsymbol{\nu}} \quad \boldsymbol{\nu} = \overline{\mathbf{1}, m}, \quad \boldsymbol{\bar{\jmath}} = \mathbf{1}, 2. \end{array} \tag{73}$$

Since the set of admissible solutions of problem (63)–(66) is larger than that of Problem S, it follows that

$$I\_S^\* \le I\_C^\*. \tag{74}$$

where *I*∗ *<sup>C</sup>* denotes the value of the optimal cyclic mode problem.

One of the problems in designing cyclic processes consists of distinguishing a class of problems for which inequality (74) turns into an equality, i.e., the cyclic extension is equivalent to the static problem. An important role in solving this problem is played by the Lagrangian function of Problem S,

$$R\_{\mathbb{S}} = f\_0(\mathbf{x}, \boldsymbol{\mu}, a) + \sum\_{\boldsymbol{\nu}} \lambda\_{\boldsymbol{\nu}} f\_{\boldsymbol{\nu}}(\mathbf{x}, \boldsymbol{\mu}, a) + \sum\_{\boldsymbol{j}} \xi\_{\boldsymbol{j}} \boldsymbol{\varphi}\_{\boldsymbol{j}}(\mathbf{x}, \boldsymbol{\mu}, a) \tag{75}$$

To determine whether a cyclic process is equivalent to a static one or efficient without solving problem (63)–(66), we form averaged problems, which are in turn extensions for Problem S or C or for both. Comparing the values of these problems with *I*∗ *<sup>C</sup>* , we find conditions for the equivalence of a cyclic extension.

1. *An upper bound for I*∗ *<sup>C</sup> and sufficient conditions for the equivalence of a cyclic extension*. Let us enlarge the set of admissible solutions of Problem C by removing the differential equations (63). We obtain Problem *S*, which we call an *estimating* problem:

$$I\_{\mathcal{S}} = \overline{f\_0(\mathbf{x}, u, a)}^{\mathbf{x}, \mu} \to \max \left/ \begin{array}{c} \overline{f\_{\nu}(\mathbf{x}, u, a)}^{\mathbf{x}, \mu} = 0, \quad \overline{q\_j(\mathbf{x}, u, a)}^{\mathbf{x}, \mu} = 0 \\ \nu = \overline{1, m}, \quad u \in V\_{\mathrm{ul}}, \quad a \in V\_a, \quad j = \overline{1, r}. \end{array} \right. \tag{76}$$

Clearly,

$$I\_{\mathcal{S}}^\* \ge I\_{\mathcal{C}'}^\* \tag{77}$$

and Problem *S* is an averaged extension of Problem S with the variables *x* and *u* and the parameters *a*. The roles of the variables *x* and *u* in the conditions of Problem *S* are similar, and we unite these variables and denote them by *y* = (*x*, *u*). In shorthand notation, this problem has the form

$$I\_{\mathbb{S}} = \overline{f\_0(y, a)}^y \to \max \left/ \overline{f\_\nu(y, a)}^y = 0, \quad \overline{\overline{g\_j(y, a)}}^y = 0, \quad \nu = \overline{1, m}, \quad j = \overline{1, r}. \tag{78}$$

The value of problem (78) as an extension of the optimal static mode problem can be expressed in terms of the function *RS* as

$$I\_{\mathbb{S}}^{\*} = \inf\_{\lambda, \mathbb{S}} \sup\_{y} R\_{\mathbb{S}}\left(y, a^{\*}, \lambda, \mathbb{S}\right). \tag{79}$$

For determining the vector of parameters, we have the condition

$$\left[\frac{\partial}{\partial a}\overline{R\_S(y,a,\lambda,\xi)}^y\right]\_{a=a^\*} = 0.\tag{80}$$

If *a*∗ lies inside *Va*, then condition (80) reduces to the condition of stationarity of *RS* with respect to *a*.

If the value *I*∗ <sup>5</sup> given by (79) equals *I*<sup>∗</sup> <sup>S</sup> (i.e., Problem S has a unique base solution), then inequalities (74) and (77) imply *I*∗ <sup>C</sup> = *I*<sup>∗</sup> <sup>S</sup> ; i.e., the static mode cannot be improved by passing to a cyclic mode. If *I*∗ *<sup>S</sup>* <sup>&</sup>gt; *<sup>I</sup>*<sup>∗</sup> <sup>S</sup> , then the difference Δ<sup>S</sup> between these values gives an upper bound for the possible gain from the passage to a cyclic mode.

2. *A lower bound for I*∗ <sup>C</sup>. *Quasi-static and sliding modes*. Consider the case when *x*(*t*) and *u*(*t*) vary so that the time derivatives of *x*(*t*) can be neglected. Then the relations between *x* and *u* are given, as in the static case, by *f*(*x*(*t*), *u*(*t*), *a*) = 0 for all *t*. The corresponding modes are said to be quasi-static. The problem of an optimal choice of *x*(*t*) and *u*(*t*) under the quasi-static conditions (Problem QS) has the form

$$I\_{\rm QS} = \frac{1}{\pi} \int\_0^\pi f\_0(\mathbf{x}, u, a) dt \to \max \, / f(\mathbf{x}, u, a) = 0, \quad \int\_0^\pi \varrho(\mathbf{x}, u, a) dt = 0, \quad u \in V\_{\rm H}, \quad a \in V\_{\rm H}.$$

or, in shorthand notation,

$$I\_{\rm QS} = \overline{f\_0(y, a)}^y \to \max \sqrt{\rho(y, a)}^y = 0, \quad y \in V\_{y\prime} \quad a \in V\_a. \tag{81}$$

Here *y* = (*x*, *u*), and the set *Vy* is determined by the conditions *u* ∈ *Vu*, *a* ∈ *Va*, and *f*(*x*, *u*, *a*) = 0. Since any solution of Problem QS is admissible for Problem C, it follows that

$$I\_{\rm QS}^\* \le I\_{\rm C}^\*. \tag{82}$$

At the same time, the value *I*∗ *QS* of Problem QS, being the value of an averaged problem, is given by the expression

$$I\_{\rm QS}^\* = \inf\_{\xi} \left\{ \sup\_{y} \left[ f\_0 \left( y, a^\* \right) + \xi \varphi \left( y, a^\* \right) \right] \, \Big/ f \left( y, a^\* \right) = 0, u \in V\_u \right\}. \tag{83}$$

Here *a*∗ is the optimal value of *a* subject to the constraint

$$\frac{\partial}{\partial a} \left[ \overline{f\_0(y, a)}^y + \left( \overline{\zeta}, \overline{\varrho(y, a)}^y \right) + \sum\_{i=0}^r \lambda^i f\left( y^i, a \right) \right]\_{a^\*} \delta a \le 0. \tag{84}$$

in which *δa* is the set of variations allowed by the inclusion *a* ∈ *Va*.

We choose the Lagrange multipliers *λ<sup>i</sup>* in (84) so that *f*(*y<sup>i</sup>* , *a*) = 0 for any base value *y<sup>i</sup>* of the vector *y*. The number of base values of *y* is determined by the dimension *r* of the vector function *ϕ*; thus, the problem takes the form

$$\bar{f}\_0 = \sum\_{i=0}^r \gamma\_i f\_0 \left( y^i, a \right), \quad \bar{\varphi} = \sum\_{i=0}^r \gamma\_i \varphi \left( y^i, a \right), \quad \sum\_{i=0}^r \gamma\_i = 1, \quad \gamma\_i \ge 0. \tag{85}$$

Consider the case when the control vector in the steady state of the system changes with a frequency so high that the state vector *x* remains virtually constant. Such a mode is called a *sliding* steady mode. The optimization problem for such a mode is formulated as

$$I\_{\rm SL} = \overline{f\_0(b, u)^u} \to \text{sup} \, / \overline{f(b, u)}^u = 0, \quad u \in V\_u, \quad \overline{\varrho(b, u)}^u = 0, \quad b \in V\_b. \tag{86}$$

This problem is known as Problem SL. In (86), *b* denotes the vector formed by *x* and *a*. This mode is the limit case of the cyclic mode, so we have

$$I\_{\mathbb{S}\mathbb{L}}^\* \le I\_{\mathbb{C}^\*}^\*$$

Problem (86) is an averaged extension of Problem S with two types of variables; its value is given by

$$I\_{\rm SL}^{\*} = \min\_{\lambda, \mathfrak{z}\_{\rm \mathfrak{z}}^{\rm x}} \max\_{u} R\left(u, b^{\*}, \lambda, \mathfrak{z}\_{\rm \mathfrak{z}}^{\rm x}\right) = \min\_{\lambda, \mathfrak{z}\_{\rm \mathfrak{z}}^{\rm x}} \max\_{u \in V\_{\rm \mathfrak{z}}} \left[f\_{0}\left(u, b^{\*}\right) + \lambda f\left(u, b^{\*}\right) + \xi \,\upvarphi\left(u, b^{\*}\right)\right],\tag{87}$$

where *b*∗ satisfies the condition

$$\frac{\partial}{\partial b}\overline{\mathcal{R}\left(\mu,b^\*,\lambda,\mathcal{J}\right)}^{\mu}\delta b \le 0. \tag{88}$$

The number of base values of the vector function *u* in Problem SL is at most *m* + *r* + 1.

A necessary condition for the efficiency of the transition to a cyclic mode can be stated in terms of *I*∗ QS and *I*<sup>∗</sup> SL. Consider the quantity

$$I\_{\mathbb{K}} = \max\left[I\_{\mathbb{QS}'}^{\*}I\_{\mathbb{SL}}^{\*}\right].$$

If *I*<sup>K</sup> is greater than *I*∗ <sup>S</sup> , then the passage to a cyclic mode is efficient, and the difference

$$
\Delta\_{\mathbb{K}} = I\_{\mathbb{K}} - I\_{\mathbb{S}}^\*
$$

provides a lower bound for the efficiency.

## *5.2. The Frequency Criterion for the Efficiency of the Passage to a Cyclic Mode*

Suppose that an optimal static mode *x*0, *u*<sup>0</sup> in Problem S is known. As above, it is required to determine whether the cyclic extension of Problem S is efficient. In [14], a frequency criterion for the efficiency of a cyclic mode was proposed. This criterion is based on the analysis of the increment in the optimality criterion *I* as compared to its maximal static value *I*<sup>0</sup> for small harmonic oscillations of the control about *u*0.

Let *λ*<sup>0</sup> and *μ*<sup>0</sup> be the values of Lagrange multipliers *λ* and *μ* corresponding to the optimal static mode in the Lagrangian function

$$R = f\_0(\mathbf{x}, \boldsymbol{\mu}) + \sum\_{i} \lambda\_i f\_i(\mathbf{x}, \boldsymbol{\mu}) + \sum\_{j} \mu\_j \boldsymbol{\uprho}\_j(\mathbf{x}, \boldsymbol{\mu})$$

for Problem S.

In a neighborhood of the optimal static mode and the corresponding Lagrange multipliers, we calculate the first and second derivatives of the functions that determine the problem with respect to *x* and *u* (if *x* and *u* are vectors, then these derivatives are matrices):

$$A = \frac{\partial f}{\partial \mathbf{x}'}, \quad B = \frac{\partial f}{\partial \mathbf{u}'}, \quad P = \frac{\partial^2 R}{\partial \mathbf{x}^2}, \quad Q = \frac{\partial^2 R}{\partial \mathbf{x} \partial \mathbf{u}'}, \quad H = \frac{\partial^2 R}{\partial \mathbf{u}^2}, \quad K = \frac{\partial \boldsymbol{\varrho}}{\partial \mathbf{x}'}, \quad M = \frac{\partial \boldsymbol{\varrho}}{\partial \mathbf{u}}.$$

In a neighbourhood of the optimal static mode, the increment of the functional *I* under small variations *δx*(*t*) and *δu*(*t*) is given by

$$
\Delta I = \frac{1}{2T} \int\_0^T \left( \delta \mathbf{x'} P \delta \mathbf{x} + \delta \mathbf{x'} Q \delta u + \delta u' Q' \delta \mathbf{x} + \delta u' H \delta u \right) dt.
$$

The transition to a cyclic mode is efficient if there is a variation *δu* such that the quantity Δ*I* is positive under the linearized constraints (63), i.e.,

$$
\delta \dot{\mathbf{x}} = A \delta \mathbf{x} + B \delta u, \quad \delta \mathbf{x}(T) = \delta \mathbf{x}(0). \tag{89}
$$

To get rid of these constraints, we consider only harmonic variations, i.e., those of the form

$$\delta\mu(t) = \sum\_{\nu=-\infty}^{\infty} \mu\_{\nu} e^{i\nu \frac{2\pi}{T}t}.$$

Applying the Fourier transform to the linear differential constraints (89), we obtain

$$
\delta\mathfrak{x}(i\omega) = \delta\mathfrak{u}(i\omega) \frac{B}{i\omega E - A} = \delta\mathfrak{u}(i\omega)\mathcal{W}(i\omega).
$$

Here *E* is the identity matrix (*E* = 1 for scalar *x*). It is assumed that the matrix *A* has no eigenvalues with zero real part; otherwise, small deviations *δu*(*t*) may correspond to large deviations *δx*(*t*), and the linearization may be incorrect.

Let us express the quantity Δ*I* by Parseval's identity in the frequency domain, replacing *δx*(*iω*) by its expression in terms of *δu*. The increment in the criterion under harmonic oscillations of the control with frequencies that are multiples of 2*π*/*T* takes the form

$$
\Delta I = \frac{1}{2} \int\_{-\infty}^{\infty} \delta u'(-i\omega) A(\omega) \delta u(i\omega) d\omega.
$$

Here *A*(*ω*) is defined by the matrices *P*, *Q*, and *H* and the relation between *δu* and *δx*; it is easy to show that

$$A(\omega) = \mathcal{W}'(-i\omega)P\mathcal{W}(i\omega) + \mathcal{Q}'\mathcal{W}(i\omega) + \mathcal{W}'(-i\omega)Q + H\_{\omega}$$

where the prime denotes transposition.

For the scalar problem, we have

$$A(\omega) = P|\mathcal{W}(i\omega)|^2 + 2Q \operatorname{Re} \mathcal{W}(i\omega) + H$$

If the matrix *A*(*ω*) for some *ω* is such that the integrand in the expression for Δ*I* is positive for at least one vector *δu*, then the static mode can be improved and the passage to a cyclic mode is efficient.

For the scalar problem, we have

$$
\Delta I = \frac{1}{2} \int\_{-\infty}^{\infty} |\delta u(i\omega)|^2 A(\omega) d\omega\_{\omega}
$$

and the static mode improves if *A*(*ω*) is positive for some *ω*.

## *5.3. Lyapunov Problems*

For an important class of problems, the inequality (77) turns into an equality. In these problems, the functions *f*0, *f* , and *ϕ* in relations (63)–(66) depend only on *u* and *a*, so that

$$
\dot{\mathfrak{x}} = f(u, a) \tag{90}
$$

Such equations are called *Lyapunov-type equations*, and the corresponding problems are known as Lyapunov problems. If we discard equations (63), which have the form (90), in Problem C, thereby transitioning to Problem S, then we can find its solution *u*∗(*t*) , *a*∗. Substituting this solution into equation (90), we determine an optimal trajectory. Clearly, in this problem, *I*∗ *<sup>S</sup>*¯ = *I*<sup>∗</sup> *<sup>C</sup>*, *u*∗(*t*) takes at most *m* + *r* + 1 base values, and the function *x*∗(*t*) is a polygonal line with at most *m* + *r* (internal) vertices.

Problems that include, in addition to Lyapunov-type equations, equations of the form

$$\dot{\mathfrak{x}}\_{\mathcal{V}} = f\_{\mathcal{V}}(\mathfrak{u}, a) F\_{\mathcal{V}}(\mathfrak{x}\_{\mathcal{V}})$$

can also be reduced to Lyapunov problems. Indeed, such equations can be reduced to the form (90) by the change

$$\mathcal{Y}\_{\nu}\left(\mathbf{x}\_{\nu}\right) = \int \frac{d\mathbf{x}\_{\nu}}{F\_{\nu}\left(\mathbf{x}\_{\nu}\right)}\,,\tag{91}$$

so that *y*˙*<sup>ν</sup>* = *fν*(*u*, *a*). The optimal solution *yν*(*t*) is piecewise linear, and *xν*(*t*) can be found from (91) by solving the equation

$$\frac{d\mathbf{x}\_{\boldsymbol{\nu}}}{F\_{\boldsymbol{\nu}}(\mathbf{x}\_{\boldsymbol{\nu}})} = y\_{\boldsymbol{\nu}}(t)dt.$$

## **6. Average Optimization in Finite-Time Thermodynamics**

The field of finite-time thermodynamics is one of the most important examples of application of averaged optimization techniques. The reasons for this are the following:

1. Problems of optimal thermodynamic cycles.

There are a very important kind of thermodynamic systems — intermediary ones. These systems contact different subsystems (reservoirs) alternately while producing power and thus lowering the irreversibility arising from a continuous contact of the above-mentioned subsystems. The main example here is the heat engine, where the working fluid contacts two sources of different temperature.

One of the most essential problems in finite-time thermodynamics is the problem of maximum average power of heat engines, when the average rate of the heat flow from the hot source is given.

Similar problems arise also in absorption-desorption systems, where the working fluid contacts with the multi-component mixture and picks one component out from one source, releasing it at another one.

In reverse cycles, the working fluid obtains the energy from the exterior system. Upon contact with the source that loses energy or matter in the regular cycle, the working fluid enriches it with the corresponding resource.

In all of these problems, the working fluid restores its state at the beginning of every cycle. One needs to average all of the variables determining the process.

2. Relations between intensive and extensive variables are Lyapunov-type equations. Thermodynamic variables are divided into two classes: intensive (temperature, pressure, chemical potential, . . . ) and extensive (volume, internal energy, entropy, amount of substance, ...) ones. Flow rates of transport processes between subsystems depend only on intensive variables. This value determines in turn the rate of change of extensive variables. This means that equations determining the change of state of the thermodynamic system have the form [10–12]:

$$\frac{dZ\_{\dot{j}}}{dt} = F\_{\dot{j}}(u\_{i\prime}u\_{\dot{j}}).\tag{92}$$

Here *i* and *j* are indices of the contacting subsystems, *u* is the vector of intensive variables, *Z* is the vector of extensive variables. Equations of this type are called Lyapunov-type equations earlier in this paper. The right hand side of these equations does not depend on *Z*, and the increase of *Z* is determined by the average value of the function *F*. As we have shown above, one can obtain the limiting capabilities of systems characterized by Lyapunov-type equations using techniques of the averaged optimization.

## **7. Example: Averaged Optimization of a Heat Engine**

## *7.1. Maximum Average Power Output*

We will assume that there is a heat engine with constant temperature of sources *T*+ and *T*<sup>−</sup> [15]. If we denote the temperature of the source contacting with the working fluid at the moment as *Tn* and the temperature of the working fluid itself as *T*, we will obtain that the average output per cycle is

$$
\overline{p} = \overline{q(T\_{n\_\prime} T)}.\tag{93}
$$

Now we can formulate the averaged optimization problem, given that the average entropy generation within the working fluid per cycle is zero:

$$\overline{q(T\_n, T)} \to \max\_T \sqrt{\begin{cases} \frac{q(T\_n, T)}{T} = 0, \\ T\_n = (T\_+; T\_-), \quad T > 0. \end{cases}}\tag{94}$$

This is the problem in the form (6). Using the algorithm described earlier (see (42)–(44)) we find the number of base points is two. This means that the Lagrange function

$$R = q(T\_{n\prime}, T) + \lambda \frac{q(T\_{n\prime}, T)}{T} = q(T\_{n\prime}, T) \left(1 + \frac{\lambda}{T}\right) \tag{95}$$

has two maxima, so

$$\begin{aligned} T\_1 &= \underset{T}{\arg\max} \, q(T\_+, T) \left(1 + \frac{\lambda}{T}\right), \\ T\_2 &= \underset{T}{\arg\max} \, q(T\_-, T) \left(1 + \frac{\lambda}{T}\right). \end{aligned} \tag{96}$$

Both maxima are global and therefore they must be equal [16]. It means that the Lagrange multiplier is the solution of

$$q(T\_{+},T\_{1})\left(1+\frac{\lambda}{T\_{1}}\right) = q(T\_{+},T\_{2})\left(1+\frac{\lambda}{T\_{2}}\right).\tag{97}$$

When the heat transfer law is linear

$$q(T\_+,T) = \mathfrak{a}\_+(T\_+ - T), \quad q(T\_-,T) = \mathfrak{a}\_-(T\_- - T), \tag{98}$$

solution of equations (94)–(96) with notation *α* = *<sup>α</sup>*+*α*<sup>−</sup> ( <sup>√</sup>*α*++√*α*−) <sup>2</sup> leads to

$$p\_{\max} = \alpha \left(\sqrt{T\_+} - \sqrt{T\_-}\right)^2. \tag{99}$$

The relationship between entropy generation and heat flows is shown at Figure 3. It is clear from this picture that the point of maximum power output lies on the convex hull of original output curves, so it is attainable only when averaged control is used.

**Figure 3.** Relationship between entropy generation and heat flows and its convex hull. Here *q<sup>h</sup>* and *q<sup>c</sup>* are heat exchange rates upon contact with the hot and cold reservoirs, respectively, and *σh*, *σ<sup>c</sup>* are the corresponding entropy generation rates. The optimal solution is attained when *σ<sup>c</sup>* = *σ*1, *σ<sup>h</sup>* = *σ*2, *q*<sup>1</sup> = *qc*(*σ*1), *q*<sup>2</sup> = *qh*(*σ*2) and *pmax* = *q*<sup>1</sup> + *q*2.

## *7.2. Maximum Efficiency*

When the power output *p*<sup>0</sup> is given, the problem of maximum efficiency is equivalent to the problem of minimum entropy generation within the system. Using again that the average entropy generation within the working fluid is zero for a cyclic process, we obtain the problem:

$$-\sigma = \overline{\left(\frac{q(T\_{\rm n}, T)}{T\_{\rm n}}\right)} \rightarrow \max\_{\overline{T}} \text{,} \sqrt{\frac{\overline{q(T\_{\rm n}, T)}}{\left(\frac{q(T\_{\rm n}, T)}{T}\right)}} = 0,\tag{100}$$

$$T\_{\rm n} = (T\_{+}; T\_{-}), \quad T > 0.$$

One may notice that this problem allows three base points in general, because there are three averaging operations in (100). This is the case when the entropy generation as a function of *q*(*T*+, *T*) is not convex. We will not consider this case here, because this function is convex for most of heat transfer laws.

Another possibility corresponds to two base points. In this case, we have the following equations for *T*<sup>1</sup> and *T*2:

$$\begin{aligned} T\_1 &= \underset{T}{\arg\max} \left[ q(T\_+, T) \left( \frac{1}{T\_+} + \lambda + \frac{\mu}{T} \right) - \lambda p\_0 \right], \\ T\_2 &= \underset{T}{\arg\max} \left[ q(T\_-, T) \left( \frac{1}{T\_-} + \lambda + \frac{\mu}{T} \right) - \lambda p\_0 \right]. \end{aligned} \tag{101}$$

These maxima must be equal, which leads to:

$$q(T\_+,T)\left(\frac{1}{T\_+} + \lambda + \frac{\mu}{T}\right) = q(T\_-,T)\left(\frac{1}{T\_-} + \lambda + \frac{\mu}{T}\right).\tag{102}$$

The averaged constraints must also be satisfied:

$$\begin{cases} \gamma q(T\_+, T\_1) + (1 - \gamma) q(T\_-, T\_2) = p\_{0\prime} \\ \gamma \frac{q(T\_+, T\_1)}{T\_1} + (1 - \gamma) \frac{q(T\_-, T\_2)}{T\_2} = 0. \end{cases} \tag{103}$$

Equations (101)–(103) allow one to find values of *T*1, *T*2, *λ*, *μ* and *γ*.

For the linear heat transfer law we have the following value of maximum efficiency:

$$\eta\_{\max}(p) = \frac{1}{2} \left( \frac{p}{aT\_+} + \eta\_\varepsilon \right) \pm \sqrt{\frac{1}{4} \left( \frac{p}{aT\_+} + \eta\_\varepsilon \right)^2 - \frac{p}{aT\_+}}.\tag{104}$$

When *p* → 0, the value of (104) approaches *η<sup>c</sup>* (Carnot efficiency) and when *p* = *pmax* (99), we have

$$
\eta\_{\max}(p\_{\max}) = 1 - \sqrt{\frac{T\_{-}}{T\_{+}}} = 1 - \sqrt{1 - \eta\_{c}} \tag{105}
$$

which is the well-known result of Novikov [17], Chambadal [18], Curzon and Ahlborn [19]. Important results for other types of heat transfer laws and different processes are presented in [20–22].

## **8. Results**

We obtained the general necessary conditions of optimality for averaged optimization problems. These conditions can be written down using the algorithmic procedure given in the paper, which allows one to use them for problems of any structure. We showed how these techniques can be applied to the problems of finite-time thermodynamics leading to new results in the field.

**Author Contributions:** Conceptualization, A.T. and I.S.; methodology, A.T.; validation, I.S.; writing–original draft preparation, A.T.; writing–review and editing, I.S.; supervision, A.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

## **References**


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