Evolutionary Search for Globally Optimal Stable Multicycles in Complex Systems with Inventory Couplings

Abstract

This note is devoted to multiperiodically operated complex system with inventory couplings transferring waste products from some subsystems as useful components to other subsystems. The flexibility of the inventory couplings is used to force each of the subsystems with its own period and to exploit its particular dynamic properties. This enhances the performance of the complex system endowed with many recycling loops, which reduce the amount of waste products endangering the natural environment. The subsystems are characterized by generalized populations composed of the individuals (the cycles), each of them encompasses its period, its initial state, its local control, and its inventory interaction. An evolutionary optimization algorithm employing such generalized populations coordinated on the basis of the inventory interaction constraints is developed. It includes the stability requirements imposed on the cyclic control processes connected with particular subsystems. The algorithm proposed is applied to the global multiperiodic optimization of some interconnected chemical production processes.

1. Introduction

Optimization of complex production systems composed of cooperating subsystems may essentially enhance their overall performance. Multiperiodically operated complex system with the inventory couplings offers new possibilities of improving the system productivity. In such systems waste products of some subsystems can be transferred to suitably chosen cooperating subsystems utilizing them as useful process components. This concept is connected with the tendency of the rearrangement of complex industrial production systems from an open loop form, yielding many waste products, to a closed loop form recycling waste products, and having desired ecological features [1–3]. The couplings of inventory type ensure a high degree of autonomy of the subsystems. Their operation periods may be chosen independently to some extent, and adjusted to their particular dynamic properties. This leads to the choice of the multicycle defined as the set of the cycles, each of them encompasses the period, the initial state, the local control, and the inventory interaction connected with the particular subsystems. Thus a new approach of the multiperiodic control is applicable to the class of discussed systems, which generalizes the periodic control approach.

The optimal periodic control (OPC) problem has reached much attention in the literature, where various periodic optimization methods have been presented exploiting the gradient-type hill-climbing [4], the second variation algorithms [5–12], the Pontryagin maximum principle [9, 13–15], the dynamic programming principle [16], the asymptotic series approach [17], the basis functions expansion [18–20], the differential flatness approach [21], and the evolutionary search for globally optimal stable cycles [22]. The latter method allows us to incorporate into the process optimization various side constraints difficult for other methods. This concerns the integral state and control constraints, the pointwise state constraints, and the process stability requirements, which may be represented by the restriction of the modulus of the Floquet's multipliers. Assuming low values of these modulus we determine highly stable optimized cycles easy implementable with the help of simple regulating loop, for instance, of relay-type. Moreover, the evolutionary search is oriented at the finding of a globally optimal cycle guaranteeing the best degree of the steady-state performance improvement. This is achieved with the aid of a population of cycles evolving by the crossing, and the mutation operations imitating the mechanism of the natural selection [23, 24].

In the present paper the evolutionary algorithm described in [22] is generalized to complex systems with the inventory couplings. The cycle of each subsystem is encoded by the period duration, by the initial state, by the discretized control, and by the discretized inventory interaction. Such cycles are locally evolved on the basis of genetic optimization scheme, and globally coordinated with the aid of the averaged inventory interaction constraints. Because of the flexibility of the interaction constraints, the operation periods and the stability margins of the subsystems are chosen independently to the extent following from the coordination constraints. The generalized evolutionary algorithm is applied to the multiperiodic optimization of the complex chemical production systems performed in two cooperating cross-recycled chemical reactors. The issue of the paper is referred to the reactor network synthesis and the waste reduction [25–28].

2. Problem Formulation

Consider the following globally optimal multiperiodic control (GOMC) problem for systems composed of 𝑁-subsystems with the inventory couplings (IC): minimize the global objective function

𝐺((𝜏_𝑖,𝑥_𝑖,𝑢_𝑖,𝑣_𝑖)^𝑁_𝑖=1)≐^𝑁_𝑖=11𝜏_𝑖^𝜏𝑖₀𝑔_𝑖(𝑥_𝑖(𝑡),𝑢_𝑖(𝑡),𝑣_𝑖(𝑡))𝑑𝑡 (1) subject for 𝑖=1,2,…,𝑁 to the 𝜏_𝑖-periodic state equations of the subsystems

𝑥̇_𝑖(𝑡)=𝑓_𝑖(𝑥_𝑖(𝑡),𝑢_𝑖(𝑡),𝑣_𝑖(𝑡)),𝑡∈[0,𝜏_𝑖],𝑥_𝑖(0)=𝑥_𝑖(𝜏_𝑖), (2) to the inclusion constraints

𝜏_𝑖∈𝒯_𝑖,𝑥_𝑖(𝑡)∈𝑋_𝑖,𝑢_𝑖(𝑡)∈𝑈_𝑖,𝑣_𝑖(𝑡)∈𝑉_𝑖,𝑡∈[0,𝜏_𝑖], (3) to the averaged control constraints

1𝜏_𝑖^𝜏𝑖₀𝑢_𝑖(𝑡)𝑑𝑡=𝑏_𝑖, (4) to the stability constraints

|𝑠_𝑖(Φ_𝑖(𝜏_𝑖,𝑥_𝑖,𝑢_𝑖,𝑣_𝑖))|_∞≤𝛼_𝑖, (5) and to the averaged inventory interaction constraints

1𝜏_𝑖^𝜏𝑖₀𝑣_𝑖(𝑡)𝑑𝑡≤^𝑁_𝑗=11𝜏_𝑗^𝜏𝑗₀ℎ_𝑖𝑗(𝑥_𝑗(𝑡),𝑢_𝑗(𝑡),𝑣_𝑗(𝑡))𝑑𝑡, (6) where 𝜏_𝑖∈𝑅₊ is the operation period of the 𝑖th subsystem, 𝑥_𝑖(𝑡)∈𝑅^𝑛𝑖 is its state at time 𝑡, 𝑢_𝑖(𝑡)∈𝑅^𝑚𝑖 is its control at time 𝑡, 𝑣_𝑖(𝑡)∈𝑅^𝑟𝑖 is its inventory interaction at time 𝑡, and the hyperparallelepipeds 𝒯_𝑖≐[𝜏⁻_𝑖,𝜏⁺_𝑖],𝑋_𝑖≐[𝑥⁻_𝑖,𝑥⁺_𝑖],𝑈_𝑖≐[𝑢⁻_𝑖,𝑢⁺_𝑖],𝑉_𝑖≐[𝑣⁻_𝑖,𝑣⁺_𝑖] with the bounds 𝜏^±_𝑖∈𝑅₊,𝑥^±_𝑖∈𝑅^𝑛𝑖,𝑢^±_𝑖∈𝑅^𝑚𝑖,𝑣^±_𝑖∈𝑅^𝑟𝑖 determine the admissible values for the process variables of the 𝑖th subsystem, and the vector 𝑏_𝑖∈𝑅^𝑚𝑖 depicts its averaged resource availability or its averaged system flow and

𝑔_𝑖∶𝑅^𝑛𝑖×𝑅^𝑚𝑖×𝑅^𝑟𝑖→𝑅,𝑓_𝑖∶𝑅^𝑛𝑖×𝑅^𝑚𝑖×𝑅^𝑟𝑖→𝑅^𝑛𝑖,ℎ_𝑖𝑗∶𝑅^𝑛𝑗×𝑅^𝑚𝑗×𝑅^𝑟𝑗→𝑅^𝑟𝑖 (7) are continuous functions defined on the sets 𝑋_𝑖×𝑈_𝑖×𝑉_𝑖, and 𝑋_𝑗×𝑈_𝑗×𝑉_𝑗, respectively, while 𝑠_𝑖≐(𝑠_𝑖𝑗)^𝑛𝑖_𝑗=1 is the vector of the eigenvalues of the monodromy matrix Φ_𝑖(𝜏_𝑖,𝑥_𝑖,𝑢_𝑖,𝑣_𝑖) of the state equation (2), |𝑠_𝑖|_∞≐max_{𝑗=1,2,…,𝑛𝑖}|𝑠_𝑖𝑗|, and 𝛼_𝑖∈𝑅₊ is the local stability level of the 𝑖th subsystem.

The objective function (1) may represent the averaged production performance of the complex system or the averaged gain from the multiperiodic operation of the IC system, which takes into account the reduction of the utilization costs of the waste products. The 𝜏_𝑖-periodic state equation (2) mean that each of the subsystems may be operated with its own period suitable for its particular dynamic properties.

The constraints (5) guarantee the stability of periodic solutions according to the approach depicted in [22].

The constraints (6) mean that the averaged outflow of the 𝑖th inventory coupling cannot exceed its averaged inflow. If the initial inventory and the coupling capacity are large enough, then such the approximate balancing may be sufficient to ensure that the current inventory will be not in defficiency or in excess over long horizon of the exploitation of the IC system. On the other hand the flexibility of the constraints (6) guarantees such features of the subsystems dynamics as the operation period and the stability margin may be chosen for each subsystem independently to the degree determined by the inventory interaction constraints.

To develop an optimization algorithm suitable for the problem discussed we use the time scaling 𝑡∶=𝜏_𝑖𝑡 independently to each subsystem, which reduces the control horizon of the IC system to the unit interval [0,1]. We approximate the continuous time control of the 𝑖th subsystem 𝑢_𝑖(𝑡) and its inventory interaction 𝑣_𝑖(𝑡) with the help of the discrete time control 𝑢_𝑖(𝑡)=𝑢^𝑘_𝑖 and the inventory interaction 𝑣_𝑖(𝑡)=𝑣^𝑘_𝑖 for 𝑡∈[𝑘/𝐾,(𝑘+1)/𝐾), where 𝑢^𝑘_𝑖∈𝑅^𝑚𝑖 and 𝑣^𝑘_𝑖∈𝑅^𝑟𝑖(𝑘=0,1,…,𝐾−1). We set 𝑢^𝐾_𝑖≐(𝑢⁰_𝑖,𝑢¹_𝑖,…,𝑢_𝑖^𝐾−1) and 𝑣^𝐾_𝑖≐(𝑣⁰_𝑖,𝑣¹_𝑖,…,𝑣_𝑖^𝐾−1). We assume that the normalized nonlinear state equation of each subsystem

𝑥̇_𝑖(𝑡)=𝜏_𝑖𝑓_𝑖(𝑥_𝑖(𝑡),𝑢_𝑖(𝑡),𝑣_𝑖(𝑡)),𝑡∈[0,1] (8) has the uniquely determined solution 𝑥_𝑖(𝑡,𝜏_𝑖,𝑥_𝑖(0),𝑢_𝑖(𝑡),𝑣_𝑖(𝑡)) for every optimization argument (𝜏_𝑖,𝑥_𝑖(0),𝑢_𝑖(𝑡),𝑣_𝑖(𝑡)) satisfying the constraints (3).

We convert by this way the GOMC problem to the following normalized and discretized form: minimize the global objective function

𝐺((𝑧_𝑖)^𝑁_𝑖=1)≐^𝑁_𝑖=1¹₀𝑔_𝑖(𝑥_𝑖(𝑡,𝑧_𝑖),𝑢_𝑖(𝑡),𝑣_𝑖(𝑡))𝑑𝑡 (9) subject for 𝑖=1,2,…,𝑁 to the process periodicity constraint

𝑥_𝑖(0)−𝑥_𝑖(1,𝑧_𝑖)=0, (10) to the inclusion constraints

𝜏_𝑖∈𝒯_𝑖,𝑥_𝑖(0)∈𝑋_𝑖,𝑢^𝑘_𝑖∈𝑈_𝑖,𝑣^𝑘_𝑖∈𝑉_𝑖(𝑘=0,1,…,𝐾−1), (11) to the averaged control constraint

1𝐾^𝐾−1_𝑘=0𝑢^𝑘_𝑖=𝑏_𝑖, (12) to the pointwise state constraints

𝑥_𝑖(𝑡)∈𝑋_𝑖,𝑡∈[0,1], (13) to the stability constraint

|𝑠_𝑖(Φ_𝑖(1,𝑧_𝑖))|_∞≤𝛼_𝑖, (14) and to the averaged inventory interaction constraint

1𝐾^𝐾−1_𝑘=0𝑣^𝑘_𝑖≤^𝑁_𝑗=1¹₀ℎ_𝑖𝑗(𝑥_𝑗(𝑡,𝑧_𝑗),𝑢_𝑗(𝑡),𝑣_𝑗(𝑡))𝑑𝑡, (15) where 𝑧_𝑖≐(𝜏_𝑖,𝑥_𝑖(0),𝑢^𝐾_𝑖,𝑣^𝐾_𝑖)∈𝑅^𝑀𝑖(𝑀_𝑖≐1+𝑛_𝑖+(𝑚_𝑖+𝑟_𝑖)𝐾) is the discrete representation of a controlled cycle of the 𝑖th subsystem with the generalized control encompassing its local control, and its inventory interaction, while 𝑥_𝑖(𝑡,𝑧_𝑖) is the solution of (8) for the given 𝑧_𝑖.

Setting 𝜏_𝑖=𝜏(𝑖=1,2,…,𝑁) we reduce the above problem to the globally optimal periodic control (GOPC) problem. The choice of equal operation periods for the subsystems may facilitate the balancing of the inventory interactions. However, it does not exploit fully the dynamic properties of the subsystems possibly having the state equations of various kinds (with polynomial, rational or exponential dependence on the state and control variables). Fixing all the process variables we further reduce the discussed problem to the globally optimal steady-state (GOSS) control problem equivalent to the minimization of the global steady-state objective function

𝐺((𝑧_𝑖)^𝑁_𝑖=1)≐^𝑁_𝑖=1𝑔_𝑖(𝑥_𝑖,𝑢_𝑖,𝑣_𝑖) (16) subject for 𝑖=1,2,…,𝑁 to the steady-state constraints

0=𝑓_𝑖(𝑥_𝑖,𝑢_𝑖,𝑣_𝑖),𝑥_𝑖∈𝑋_𝑖,𝑢_𝑖∈𝑈_𝑖,𝑣_𝑖∈𝑉_𝑖,𝑢_𝑖=𝑏_𝑖,|𝑠_𝑖(𝑒^{𝑓𝑖𝑖,𝑥(𝑥𝑖,𝑢𝑖,𝑣𝑖)})|_∞≤𝛼_𝑖,𝑣_𝑖≤^𝑁_𝑗=1𝑔_𝑗(𝑥_𝑗,𝑢_𝑗,𝑣_𝑗), (17) where 𝑧_𝑖≐(𝑥_𝑖,𝑢_𝑖,𝑣_𝑖)^𝑁_𝑖=1 is the steady-state control process of the 𝑖th subsystem.

We are aimed at the comparison of the three nested kinds of the complex process operation, namely, the static one, the periodic one, and the multiperiodic one.

3. Evolutionary Algorithm for the Global Multiperiodic Optimization

The way of coding any optimization problems has a main impact on the effectiveness of the evolutionary algorithm. The wrong choice of individuals can cause that the evolutionary algorithm can stack in local solutions or time of searching can increase to not acceptable value. Owing to complexity of the GOMC problem we think that the best way of coding problem (9)–(15) is to use vectors as individuals with real numbers genes. We propose to represent an individual by the vector 𝑧≐(𝑧_𝜇)𝑁_𝜇=1, where 𝑁=∑^𝑁_𝑖=1𝑀_𝑖,  𝑧_𝑖≐𝜏_𝑖(𝑖=1,2,…,𝑁) are operation periods of all subsystems, 𝑧_{𝑗+𝑁+∑}^𝑖−1_{𝑙=1𝑛𝑙}≐𝑥_𝑖𝑗(0)(𝑗=1,2,…,𝑛_𝑖;𝑖=1,2,…,𝑁) are coordinates of initial states of all subsystems, 𝑧_{𝑘+1+(𝑗−1)𝐾+𝐾∑}^𝑖−1_{𝑙=1𝑚𝑙+𝑁+∑}^𝑁_{𝑙=1𝑛𝑙}≐𝑢^𝑘_𝑖𝑗(𝑘=0,1,…,𝐾−1;𝑗=1,2,…,𝑚_𝑖;𝑖=1,2,…,𝑁) are discrete-time control coordinates of all subsystems, 𝑧_{𝑘+1+(𝑗−1)𝐾+𝐾∑}^𝑖−1_{𝑙=1𝑟𝑙+𝐾∑}^𝑁_{𝑙=1𝑚𝑙+𝑁+∑}^𝑁_{𝑙=1𝑛𝑙}≐𝑣^𝑘_𝑖𝑗(𝑘=0,1,…,𝐾−1;𝑗=1,2,…,𝑟_𝑖;𝑖=1,2,…,𝑁) are discrete-time inventory interactions of all subsystems. Values of the individual's genes are bounded by the set 𝑍≐[𝑧⁻,𝑧⁺](𝑧^±∈𝑅𝑁), which follows from a set of inclusion constrains 𝒯×𝑋×𝑈×𝑉.

The form of the individual 𝑧 allows us to use the uniform crossing operator, which creates two new individuals 𝑧⁺¹=(𝑧⁺¹_𝜇)𝑁_𝜇=1,𝑧⁺²=(𝑧⁺²_𝜇)𝑁_𝜇=1 according to the following rule:

𝑧⁺¹_𝜇=𝑧¹_𝜇𝑧if𝑎≤0.5,_2𝜇if𝑎>0.5,₊𝑧²_𝜇=𝑧²_𝜇if𝑧⁺¹_𝜇=𝑧¹_𝜇,𝑧¹_𝜇otherwise, (18)(19) where 𝑎∈[0,1] is randomly chosen under the uniform distribution, and 𝜇=1,…,𝑁.

A nonuniform mutation operator can be used together with a uniform crossing operator. It creates a new gene 𝑧⁺_𝜇(𝜇∈{1,…,𝑁}) of the individual according to the following relationship:

𝑧⁺_𝜇=𝑧_𝜇+Δ(𝑡,̃𝑧⁺_𝜇−𝑧_𝜇)if𝑎₁𝑧≤0.5,_𝜇−Δ(𝑡,̃𝑧_𝜇−𝑧⁻_𝜇)if𝑎₁>0.5, (20) where Δ(𝑡,̃𝛿)=𝑎₂𝛿(1−𝑡/̃𝑇)^𝜎, 𝑡̃ is a number of current generation, 𝑇 is maximum number of generations, 𝑎_𝑖∈[0,1](𝑖=1,2) is chosen randomly under the uniform distribution, 𝜎>0 is the tuning parameter. The gene 𝑧_𝜇 to mutation is chosen randomly under the uniform distribution.

Crossing or mutation operators can give a new individual, which will not fulfil the periodic constraint (10), the averaged control constraint (12), and the averaged inventory interaction constraint (15).

The periodic constraint (10) can be preserved with the help of the Newton method. Using the Newton method we can calculate new initial states according to the following relationship:

𝑥⁺_𝑖(0)=𝑥_𝑖(0)−(𝐼−Φ_𝑖(1,𝑧_𝑖))⁻¹(𝑥_𝑖(0)−𝑥_𝑖(1,𝑧_𝑖))(𝑖=1,2,…,𝑁). (21) The averaged constraints (12) and (15) can be preserved with the help of the reconstruction algorithm. But before we show it we have to define the auxiliary algorithm, which creates the index set.

Algorithm 1 (the auxiliary algorithm). One has the following steps.Step 1. Create set 𝐷∈{𝑑₁,𝑑₂,…,𝑑_𝐾}, where 𝑑_𝑖=0(𝑖=1,2,…,𝐾).Step 2. Create function 𝑤(𝜒^𝑘)(𝑘=0,1,…,𝐾−1), which returns values from the range 1,2,…,𝐾 depending on the value of 𝜒^𝑘. If (1/𝐾)∑_{𝑘∈{0,1,…,𝐾−1}}𝜒^𝑘≥𝛾 the function 𝑤(𝜒^𝑘) returns the position of element 𝜒^𝑘 in the set {𝜒^𝑘}^𝐾−1_𝑘=0 sorted in nonincreasing order. Otherwise the function 𝑤(𝜒^𝑘) returns the position of element 𝜒^𝑘 in the set {𝜒^𝑘}^𝐾−1_𝑘=0 sorted in nondecreasing order.Step 3. Calculate probability 𝑝_𝑘=𝑤(𝜒^𝑘)/(∑^𝐾_𝑖=1𝑖)(𝑘=0,1,…,𝐾−1).Step 4. Calculate distribution 𝑞_𝑘=∑^𝑘_𝑗=0𝑝_𝑗(𝑘=0,1,…,𝐾−1).Step 5. Draw 𝑎∈[0,1] under the uniform distribution and next find 𝑘∈0,…,𝐾−1 for which the following condition 𝑞_𝑘−1<𝑎≤𝑞_𝑘 is fulfilled. Set 𝑑₁=𝑘.Step 6. Set 𝑙=2,𝑖=𝑑₁+1,𝑗=𝑑₁−1.Step 7. If 𝑖≥𝐾 then set 𝑖=0. If 𝑗<0 then set 𝑗=𝐾−1.Step 8. Draw 𝑎∈[0,1] under the uniform distribution. If 𝑎<0.5 go to Step 9. Otherwise go to Step 11.Step 9. If 𝑙>𝐾 then stop the algorithm. Otherwise set 𝑑_𝑙=𝑖,𝑙=𝑙+1 and 𝑖=𝑖+1.Step 10. If 𝑙>𝐾 then stop the algorithm. Otherwise set 𝑑_𝑙=𝑗,𝑙=𝑙+1,𝑗=𝑗−1 and go to Step 7.Step 11. If 𝑙>𝐾 then stop the algorithm. Otherwise set 𝑑_𝑙=𝑗,𝑙=𝑙+1,𝑗=𝑗−1.Step 12. If 𝑙>𝐾 then stop the algorithm. Otherwise set 𝑑_𝑙=𝑖,𝑙=𝑙+1 and 𝑖=𝑖+1 and go to Step 7.

Having Algorithm 1 we can define the reconstruction algorithm for an averaged control and inventory interaction constraints.

Algorithm 2 (reconstruction of control and inventory constraints). One has the following steps.Step 1 (the initialization). Select a control variable 𝑢_𝑖𝑗=(𝑢^𝑘_𝑖𝑗)^𝐾−1_𝑘=0(𝑖=1,2,…,𝑁;𝑗=1,2,…,𝑚_𝑖) or an inventory interaction variable 𝑣_𝑖𝑗=(𝑣^𝑘_𝑖𝑗)^𝐾−1_𝑘=0(𝑖=1,2,…,𝑁;𝑗=1,2,…,𝑟_𝑖) for the reconstruction (a control or inventory interaction which violates constraint (12) or (15), resp.).Step 2 (the creation of the index set). If you chose the control variable to reconstruction, set 𝜒=𝑢_𝑖𝑗 and 𝛾=𝑏_𝑖𝑗. If you chose the inventory interaction to reconstruction, set 𝜒=𝑣_𝑖𝑗 and 𝛾=𝑏̃_𝑖𝑗, where 𝑏̃_𝑖=∑^𝑁_𝑘=1∫¹₀ℎ_𝑖𝑘(𝑥_𝑘(𝑡,𝑧_𝑘),𝑢_𝑘(𝑡),𝑣_𝑘(𝑡))𝑑𝑡(𝑖=1,2,…,𝑁). Next, create the index set 𝐷∈{𝑑₁,𝑑₂,…,𝑑_𝐾} using Algorithm 1 and set 𝑙=1.Step 3 (the modification of a control or inventory interaction coordinate). If you chose the control variable, modify an element 𝑢^𝑑𝑙_𝑖𝑗 by the formula 𝑢^𝑑𝑙_𝑖𝑗=Sat(𝐾𝑏_𝑖𝑗−_{𝑘∈𝐷⧵{𝑑𝑙}}𝑢^𝑘_𝑖𝑗;𝑢⁻_𝑖𝑗,𝑢⁺_𝑖𝑗). (22) If you chose the inventory interaction variable, modify an element 𝑣^𝑑𝑙_𝑖𝑗 by the formula 𝑣^𝑑𝑙_𝑖𝑗=Sat(𝐾𝑏̃_𝑖𝑗−_{𝑘∈𝐷⧵{𝑑𝑙}}𝑣^𝑘_𝑖𝑗;𝑣⁻_𝑖𝑗,𝑣⁺_𝑖𝑗). (23) The saturation function Sat is defined as Sat(𝑦;𝑦⁻,𝑦⁺𝑦)=⁻if𝑦<𝑦⁻,𝑦if𝑦⁻≤𝑦≤𝑦⁺,𝑦⁺if𝑦⁺<𝑦. (24)Step 4 (the stop criterion). If the constraint (12) or (15), respectively, is satisfied, stop the reconstruction algorithm. Else set 𝑙=𝑙+1 and go to Step 3.

Algorithm 2 is not sufficient in order to have the averaged inventory interaction constraint (15) fulfilled. When we apply Algorithm 2 for a current individual, we receive a new vector of inventory interaction. This new vector causes that the right side of inequality (15) has also a new value. Thus, when we apply algorithm we can receive individual, which still violates the constraint (15). Therefore, apart Algorithm 2 we need another way to keep the inventory interaction constraint fulfilled. We propose to use a penalty term to take into account the mentioned constraint (15) together with state (13) and stability (14) constraints.

We introduce the following extended performance index used at the selection step [23, 24]:

𝐺(𝑧)=̇𝐺(𝑧)+𝜌₁^𝑁_𝑖=1𝜂₁(𝑥_𝑖(𝑡))+𝜌₂^𝑁_𝑖=1𝜂₂(𝑠_𝑖)+𝜌₃𝜂₃, (25) where

𝜂₁(𝑥_𝑖(𝑡))=̇max_𝑡∈[0,1)|𝑥⁺_𝑖−𝑥_𝑖(𝑡)|_∞,for𝑥_𝑖(𝑡)>𝑥⁺_𝑖,|𝑥⁻_𝑖−𝑥_𝑖(𝑡)|_∞,for𝑥_𝑖(𝑡)<𝑥⁻_𝑖,0,for𝑥_𝑖𝜂(𝑡)∈𝑋,₂(𝑠_𝑖)=̇|𝑠_𝑖(Φ_𝑖(1,𝑧_𝑖))|_∞if|𝑠_𝑖(Φ_𝑖(1,𝑧_𝑖))|_∞>𝛼_𝑖,𝜂0,otherwise,₃=̇^𝑁_𝑖=1𝑒^𝑇_𝑖(1𝐾^𝐾−1_𝑘=0𝑣^𝑘_𝑖−𝑏̃_𝑖), (26)𝜌_𝑖(𝑖=1,2,3) is a positive constant satisfying the inequality 𝜌_𝑖𝜂_𝑖>𝐺(𝑧) for 𝜂_𝑖≠0 and all cycles 𝑧 in the current population, 𝑏̃_𝑖=∑^𝑁_𝑘=1∫¹₀ℎ_𝑖𝑘(𝑥_𝑘(𝑡,𝑧_𝑘),𝑢_𝑘(𝑡),𝑣_𝑘(𝑡))𝑑𝑡(𝑖=1,2,…,𝑁),𝑒_𝑖∈𝑅^𝑟𝑖 and 𝑒_𝑖𝑗=0(𝑖=1,2,…,𝑁;𝑗=1,2,…,𝑟_𝑖)if(1/𝐾)∑^𝐾−1_𝑘=0𝑣^𝑘_𝑖𝑗≤𝑏̃_𝑖𝑗 otherwise 𝑒_𝑖𝑗=1.

Having in mind the above remarks, we propose the following.

Algorithm 3 (the evolutionary algorithm for the GOMC problem). One has the following steps.Step 1 (the initialization). Assume initial data: 𝑃 the number of individuals (cycles) in the population, the parameter 𝛽<𝑁 determining the number of individuals, which take part in competition at the selection step, 𝑇 the maximum number of generations, 𝑡_𝑠̃ and 𝜀≥0 parameters of the stop criterion, the tuning parameter 𝜎>0 (20), and the mutation probability 𝑝_𝑚. Sample randomly 𝑃 individuals under uniform distribution from the set (11) of the range constraints. Set 𝑡̃=1 for the index of the generation.Step 2 (the control reconstruction). Exploiting Algorithm 2, and the Newton method (21) reconstruct the population to 𝑃 individuals satisfying the constraints (10) and (12). If the Newton method is divergent create randomly a new individual according to Step 1, and repeat Step 2.Step 3 (the inventory interaction reconstruction). Exploiting Algorithm 2, reconstruct the population to 𝑃 individuals satisfying the constraints (15). Next, use the Newton method. If the Newton method is divergent, create randomly a new individual according to Step 1, and repeat Step 2.Step 4 (the evaluation). Evaluate the performance index 𝐺(𝑧^𝜈) for all individuals in the population (𝜈=1,2,…,𝑃), and find an individual 𝑧^min satisfying constraints (10)–(15) with the minimal value of the function 𝐺.Step 5 (the stop criterion). Let 𝑧^min(𝑡)̃ be an individual with the minimal value of the function 𝐺 at the generation 𝑡.̃ If 𝑡̃=𝑇 or |𝐺(𝑧^min(𝑡̃−𝑡_𝑠̃))−𝐺(𝑧^min(𝑡)̃)|≤𝜀, take 𝑧^min(𝑡)̃ as an optimal solution, and stop the computations.Step 6 (the selection). Repeat the following activities until the current population is empty: (i)draw from the current population 𝛽 individuals and (using the formula (25)) compute for them the modified performance index 𝐺(𝑧^𝜈)(ii)take from 𝛽 individuals to a new (𝑡̃+1) population one individual with the smallest value of the 𝐺(𝑧^𝜈)(iii)remove from the current population 𝛽 individuals. Set 𝑡̃∶=𝑡̃+1.Step 7 (the reproduction). Applying the crossing operator (18)-(19) to individuals from Step 6, retrieve the population composed of 𝑃 individuals.Step 8 (the mutation). For each cycle 𝑧^𝜈(𝜈=1,2,…,𝑃) sample randomly a number 𝑝∈[0,1] under the uniform probability distribution. If 𝑝≤𝑝_𝑚 use the mutation operator (20) to create a mutated cycle 𝑧^+𝜈. Next go to Step 2.

Theory says that the probability of finding a globally optimal solution is close to 1 when the number of generations and individuals is infinity [23, 24]. But this condition can not be fulfilled in practice. Thus we need a mechanism, which can stop an evolutionary algorithm in reasonable time. We decided to stop the algorithm after 𝑇 generations or if the performance index is not improved more than 𝜀 within 𝑡_𝑠̃ generations.

4. Examples

Example 1. Let the parallel chemical reactions 𝐴_𝑖→𝐵_𝑖,𝐴_𝑖→𝐶_𝑖 take place in two continuously stirred tank reactors with the inventory interactions (see Figure 1), where 𝐴_𝑖 is the raw material of the 𝑖th reactor, 𝐵_𝑖 is its desired product, and 𝐶_𝑖 is its waste product (𝑖=1,2). Assume that the 𝑖th reactor is 𝜏_𝑖-periodically operated, and denote by 𝑥_𝑖1(𝑡),𝑥_𝑖2(𝑡),𝑥_𝑖3(𝑡) its concentrations of 𝐴_𝑖,𝐵_𝑖,𝐶_𝑖, respectively, by 𝑢_𝑖1(𝑡) its input concentration of 𝐴_𝑖, by 𝑢_𝑖2(𝑡) its flow rate, by 𝑣_𝑖(𝑡) its inventory interaction transferring the waste product of the cooperating subsystem as the supplement of the raw material of the 𝑖th subsystem. Consider the following GOMC problem for the discussed system: minimize the global objective function 𝐺((𝜏_𝑖,𝑥_𝑖,𝑢_𝑖,𝑣_𝑖)²_𝑖=1)≐−²_𝑖=11𝜏_𝑖^𝜏𝑖₀𝑐_𝑖𝑢_𝑖2(𝑡)𝑥_𝑖2(𝑡)𝑑𝑡 (27) subject for 𝑖=1,2 to the constraints 𝑥̇_𝑖1(𝑡)=𝑢_𝑖2(𝑡)(𝑢_𝑖1(𝑡)+𝑣_𝑖(𝑡)−𝑥_𝑖1(𝑡))−𝜅_𝑖1𝑥^𝑝𝑖1_𝑖1(𝑡)−𝜅_𝑖2𝑥^𝑝𝑖2_𝑖1𝑥̇(𝑡),_𝑖2(𝑡)=−𝑢_𝑖2(𝑡)𝑥_𝑖2(𝑡)+𝜅_𝑖1𝑥^𝑝𝑖1_𝑖1𝑥̇(𝑡),_𝑖3(𝑡)=−𝑢_𝑖2(𝑡)𝑥_𝑖3(𝑡)+𝜅_𝑖2𝑥^𝑝𝑖2_𝑖1𝜏(𝑡),_𝑖∈[0.1,20],0≤𝑥_𝑖𝑘(𝑡)(𝑘=1,2,3),0≤𝑢_𝑖𝑗(𝑡)≤2,0≤𝑣_𝑖(𝑡)≤2(𝑗=1,2),𝑡∈[0,𝜏_𝑖1],𝜏_𝑖^𝜏𝑗₀𝑢_𝑖𝑗(𝑡)𝑑𝑡=1(𝑗=1,2),|𝑠_𝑖(Φ_𝑖(𝜏_𝑖,𝑥_𝑖,𝑢_𝑖,𝑣_𝑖))|_∞1≤0.5,𝜏₁^𝜏1₀𝑣₁1(𝑡)𝑑𝑡≤𝜏₂^𝜏2₀𝑥₂₃1(𝑡)𝑑𝑡,𝜏₂^𝜏2₀𝑣₂1(𝑡)𝑑𝑡≤𝜏₁^𝜏1₀𝑥₁₃(𝑡)𝑑𝑡, (28) where 𝑥_𝑖𝑗(0)=𝑥_𝑖𝑗(𝜏_𝑖)(𝑗=1,2,3),𝑥_𝑖(𝑡)≐(𝑥_𝑖1(𝑡),𝑥_𝑖2(𝑡),𝑥_𝑖3(𝑡))^𝑇 is the state of the 𝑖th subsystem, and 𝑢_𝑖(𝑡)≐(𝑢_𝑖1(𝑡),𝑢_𝑖2(𝑡))^𝑇 is its control. The reactions obey the power law with the exponents 𝑝_𝑖𝑗. The optimization goal is equivalent to the maximization of the summary gain of the cooperating systems for the gain coefficients 𝑐_𝑖 connected with the 𝑖th useful product 𝐵_𝑖(𝑖=1,2). We compare the optimal static control process with periodic, and multiperiodic control processes putting the special emphasis on the choice of the operation periods for the particular subsystems.
The evolutionary algorithm is time consuming. Therefore one should first run 𝜋-test [5, 8, 29] to check if there exists a periodic or multiperiodic solution, which can improve the efficiency of the process in comparison to steady-state solution. Assuming the following subsystem parameters 𝑐_𝑖=1,𝑛_𝑖=3,𝑚_𝑖=2,𝑝_𝑖1=2,𝑝_𝑖2=1(𝑖=1,2),𝜅₁₁=2,𝜅₁₂=0.5,𝜅₂₁=10,𝜅₂₂=3 we obtained the 𝜋-curves shown in Figure 2. The presented 𝜋-curves were calculated for optimal steady-state solution 𝑢₁=𝑢₂=(1,1)^𝑇,𝑣₁=0.64,𝑣₂=0.30,𝑥₁=(0.62,0.73,0.30)^𝑇,𝑥₂=(0.21,0.45,0.64)^𝑇 which gives the performance index 𝐺=−1.18. The form of the 𝜋-curves shows that there exists a periodic and multiperiodic control ensuring a performance index better than in case of using an optimal steady-state control. The solutions obtained with the help of evolutionary Algorithm 3, which was implemented in Mathematica system [30], confirm the results of the 𝜋 test. Algorithm 3 after 565 generations calculated the suboptimal periodic solution which gave the performance index 𝐺=−3.77. This performance index corresponds to a suboptimal period 𝜏=7.4. Algorithm 3 had to generate 837 generations in order to find a suboptimal multiperiodic solution shown in Figures 3, 4, 5, 6, and 7. The suboptimal multiperiodic controls and inventory interactions shown in Figures 3 and 4 together with 𝜏₁=8.1 and 𝜏₂=13.7 ensured the performance index 𝐺=−3.88. One can see that Algorithm 3 found a suboptimal multiperiodic solution, which improves the optimal steady-state performance index about 228% and a suboptimal periodic performance index about 3%. The presented results were obtained for the following parameters of the evolutionary Algorithm 3: 𝑃=100,𝛽=2,𝑇=2500,𝑡_𝑠̃=50,𝜖=10⁻⁶,𝜎=5,𝑝_𝑚=0.25.
Since the nonzero conditions for the product concentrations occur in the solutions obtained (see Figure 5), the start-up phase of the chemical reactors should be used to carry out the initial state with the non-zero condition on the reactant but with the zero conditions on the product concentrations to the desired state with non-zero conditions for all the process components. Next the (multi) periodic control is implemented. The start-up phase may be neglected in the system performance evaluation because its duration is small as compared with long horizon exploitation of the (multi) periodic control. The start-up phase can be implemented with the help of the control determined from the standard (minimum-time) two-point boundary value problem, for which the Pontryagin maximum principle is applicable.

Figure 1 

Diagram of the parallel chemical reactions from Example 1.

Figure 2 

The results of the 𝜋-test for Example 1.

Figure 3 

The suboptimal multiperiodic control 𝑢𝑖(𝑖=1,2) for Example 1.

Figure 4 

The suboptimal multiperiodic inventory interaction 𝑣𝑖(𝑖=1,2) for Example 1.

Figure 5 

The suboptimal multiperiodic state 𝑥1 for Example 1.

Figure 6 

The suboptimal multiperiodic state 𝑥2 for Example 1.

Figure 7 

The suboptimal multiperiodic cycle for Example 1.

Example 2. Let us consider the same GOMC problem as in Example 1 but with different interpretation of the inventory interactions. We assume that the inventory interaction 𝑣₁(𝑡) of the first subsystem transfers the waste product of the cooperating subsystem as the catalyst for the chemical reaction 𝐴₁→𝐵₁,𝐴₁→𝐶₁, while the inventory interaction 𝑣₂(𝑡) transfers the waste product of the cooperating subsystem as the supplement of its raw material. This means that the state equations of the subsystems take the form 𝑥̇₁₁(𝑡)=𝑢₁₂(𝑡)(𝑢₁₁(𝑡)−𝑥₁₁(𝑡))−𝜅₁₁𝑣^𝑞11₁(𝑡)𝑥^𝑝11₁₁(𝑡)−𝜅₁₂𝑥^𝑝12₁₁𝑥̇(𝑡),₁₂(𝑡)=−𝑢₁₂(𝑡)𝑥₁₂(𝑡)+𝜅₁₁𝑣^𝑞11₁(𝑡)𝑥^𝑝11₁₁𝑥̇(𝑡),₁₃(𝑡)=−𝑢₁₂(𝑡)𝑥₁₃(𝑡)+𝜅₁₂𝑥^𝑝12₁₁𝑥̇(𝑡),₂₁(𝑡)=𝑢₂₂(𝑡)(𝑢₂₁(𝑡)+𝑣₂(𝑡)−𝑥₂₁(𝑡))−𝜅₂₁𝑥^𝑝21₂₁(𝑡)−𝜅₂₂𝑥^𝑝22₂₁𝑥̇(𝑡),₂₂(𝑡)=−𝑢₂₂(𝑡)𝑥₂₂(𝑡)+𝜅₂₁𝑥^𝑝21₂₁𝑥̇(𝑡),₂₃(𝑡)=−𝑢₂₂(𝑡)𝑥₂₃(𝑡)+𝜅₂₂𝑥^𝑝22₂₁(𝑡), (29) where 𝑥_𝑖𝑗(0)=𝑥_𝑖𝑗(𝜏_𝑖)(𝑗=1,2,3). The reactions obey the power law with the exponents 𝑝_𝑖𝑗(𝑖=1,2;𝑗=1,2).
The results presented below were obtained for the following subsystem parameters: 𝑛_𝑖=3,𝑚_𝑖=2,𝑝_𝑖1=2,𝑝_𝑖2=1(𝑖=1,2),𝑞₁₁=2,𝜅₁₁=4,𝜅₁₂=0.8,𝜅₂₁=40,𝜅₂₂=2. Just like in the previous example the optimal static control process was compared with periodic, and multiperiodic control processes. The optimal steady-state solution 𝑢₁=𝑢₂=(1,1)^𝑇,𝑣₁=0.30,𝑣₂=0.40,𝑥₁=(0.50,0.09,0.40)^𝑇,𝑥₂=(0.15,0.94,0.30)^𝑇 gives the performance index 𝐺=−1.03. The 𝜋-test results obtained for the optimal steady-state solution (see Figure 8) show that periodic and multiperiodic processes can improve the efficiency of reactions. The results of 𝜋-test confirm results obtained by means of the evolutionary Algorithm 3. Algorithm 3 found a suboptimal periodic solution (with a period 𝜏=17.3) ensuring 𝐺=−3.62, which improves a steady-state solution about 251%. Of course Algorithm 3 also found a suboptimal multiperiodic process (see Figures 9, 10, 11, 12, 13, 14, and 15) which gave the best performance index. The controls and inventory interactions shown in Figure 9, 10, 11, and 12 and the periods 𝜏₁=19.7,𝜏₂=6.5 gave the performance index 𝐺=−3.95. The suboptimal multiperiodic process improves the optimal static process about 283% and the suboptimal periodic process about 9%. The evolutionary Algorithm 3 was run for the same parameters as in Example 1. We have to add that the suboptimal periodic solution was found after 368 generations and the suboptimal multiperiodic solution was found after 1290 generations.

Figure 8 

The results of the 𝜋-test for Example 2.

Figure 9 

The suboptimal multiperiodic control 𝑢1 for Example 2.

Figure 10 

The suboptimal multiperiodic control 𝑢2 for Example 2.

Figure 11 

The suboptimal multiperiodic inventory interaction 𝑣1 for Example 2.

Figure 12 

The suboptimal multiperiodic inventory interaction 𝑣2 for Example 2.

Figure 13 

The suboptimal multiperiodic state 𝑥1 for Example 2.

Figure 14 

The suboptimal multiperiodic state 𝑥2 for Example 2.

Figure 15 

The suboptimal multiperiodic cycle for Example 2.

5. Conclusion

Multiperiodically operated complex system with the inventory couplings has been modelled and characterized as a system with a high degree of freedom. The operation period, the initial state, the discretized control, and the discretized inventory interaction have been assumed as the local finite-dimensional optimization argument (the controlled cycle of the subsystem). The set of the cycles connected with all the subsystems has been treated as the controlled multicycle of the complex system and optimized by the evolutionary algorithm dealing with the generalized populations coordinated with the aid of the averaged inventory interaction constraints. The use of different operation periods and different stability margins of the subsystems has been taken into account. It has been shown that the multiperiodic operation of the complex cross-recycled chemical production systems may improve its global performance as compared with the static operation, and the periodic operation. This may enhance the utilization of waste products endangering the natural environment.