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

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.


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][2][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][6][7][8][9][10][11][12], the Pontryagin maximum principle [9,[13][14][15], the dynamic programming principle [16], the asymptotic series approach [17], the basis functions expansion [18][19][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, International Journal of Chemical Engineering the evolutionary search is oriented at the finding of a globally optimal cycle guaranteeing the best degree of the steadystate 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][26][27][28].

Problem Formulation
Consider the following globally optimal multiperiodic control (GOMC) problem for systems composed of Nsubsystems with the inventory couplings (IC): minimize the global objective function (1) subject for i = 1, 2, . . ., N to the τ i -periodic state equations of the subsystems to the inclusion constraints to the averaged control constraints to the stability constraints and to the averaged inventory interaction constraints where τ i ∈ R + is the operation period of the ith subsystem, x i (t) ∈ R ni is its state at time t, u i (t) ∈ R mi is its control at time t, v i (t) ∈ R ri is its inventory interaction at time t, and the hyperparallelepipeds determine the admissible values for the process variables of the ith subsystem, and the vector b i ∈ R mi depicts its averaged resource availability or its averaged system flow and are continuous functions defined on the sets X i ×U i ×V i , and 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 τ i -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 ith 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 t := τ i t 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 ith subsystem u i (t) and its inventory interaction v i (t) with the help of the discrete time control u i (t) = u k i and the inventory interaction v i (t ).We assume that the normalized nonlinear state equation of each subsystem has the uniquely determined solution x i (t, τ i , x i (0), u i (t), v i (t)) for every optimization argument (τ i , x i (0), u i (t), v i (t)) satisfying the constraints (3).
We convert by this way the GOMC problem to the following normalized and discretized form: minimize the global objective function subject for i = 1, 2, . . ., N to the process periodicity constraint to the inclusion constraints to the averaged control constraint to the pointwise state constraints to the stability constraint and to the averaged inventory interaction constraint where is the discrete representation of a controlled cycle of the ith subsystem with the generalized control encompassing its local control, and its inventory interaction, while x i (t, z i ) is the solution of (8) for the given z i .
Setting τ i = τ (i = 1, 2, . . ., N) 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 subject for i = 1, 2, . . ., N to the steady-state constraints where is the steady-state control process of the ith 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.

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 z .
Values of the individual's genes are bounded by the set Z .
The form of the individual z allows us to use the uniform crossing operator, which creates two new individuals according to the following rule: where a ∈ [0, 1] is randomly chosen under the uniform distribution, and μ = 1, . . ., N.

International Journal of Chemical Engineering
A nonuniform mutation operator can be used together with a uniform crossing operator.It creates a new gene + z μ (μ ∈ {1, . . ., N}) of the individual according to the following relationship: where ) is chosen randomly under the uniform distribution, σ > 0 is the tuning parameter.The gene z μ 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: 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 5. Draw a ∈ [0, 1] under the uniform distribution and next find k ∈ 0, . . ., K − 1 for which the following condition Step 7. If i ≥ K then set i = 0.If j < 0 then set j = K − 1.
Step 8. Draw a ∈ [0, 1] under the uniform distribution.If a < 0.5 go to Step 9. Otherwise go to Step 11.
Step 9.If l > K then stop the algorithm.Otherwise set d l = i, l = l + 1 and i = i + 1.
Step 10.If l > K then stop the algorithm.Otherwise set d l = j, l = l + 1, j = j − 1 and go to Step 7.
Step 11.If l > K then stop the algorithm.Otherwise set d l = j, l = l + 1, j = j − 1.
Step 12.If l > K then stop the algorithm.Otherwise set d l = i, l = l + 1 and i = i + 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 2 (the creation of the index set).If you chose the control variable to reconstruction, set χ = u i j and γ = b i j .If you chose the inventory interaction to reconstruction, set χ = v i j and γ = b i j , where b i = N k=1 Step 3 (the modification of a control or inventory interaction coordinate).If you chose the control variable, modify an element u dl i j by the formula If you chose the inventory interaction variable, modify an element v dl i j by the formula The saturation function Sat is defined as International Journal of Chemical Engineering 5 Step 4 (the stop criterion).If the constraint ( 12) or ( 15), respectively, is satisfied, stop the reconstruction algorithm.
Else set l = l + 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]: where ) is a positive constant satisfying the inequality ρ i η i > G(z) for η i / = 0 and all cycles z in the current population, 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: P the number of individuals (cycles) in the population, the parameter β < N determining the number of individuals, which take part in competition at the selection step, T the maximum number of generations, t s and ε ≥ 0 parameters of the stop criterion, the tuning parameter σ > 0 (20), and the mutation probability p m .Sample randomly P individuals under uniform distribution from the set (11) of the range constraints.Set t = 1 for the index of the generation.
Step 2 (the control reconstruction).Exploiting Algorithm 2, and the Newton method (21) reconstruct the population to P 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 P 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 G( z ν ) for all individuals in the population (ν = 1, 2, . . ., P), and find an individual z min satisfying constraints ( 10)- (15) with the minimal value of the function G.
Step 5 (the stop criterion).Let z min ( t) be an individual with the minimal value of the function G at the generation t.If t = T or |G( z min ( t − t s )) − G( z min ( t))| ≤ ε, take z min ( t) as an optimal solution, and stop the computations.
Step 6 (the selection).Repeat the following activities until the current population is empty: Set t := t + 1.
Step 7 (the reproduction).Applying the crossing operator ( 18)-( 19) to individuals from Step 6, retrieve the population composed of P individuals.
Step 8 (the mutation).For each cycle z ν (ν = 1, 2, . . ., P) sample randomly a number p ∈ [0, 1] under the uniform probability distribution.If p ≤ p m use the mutation operator (20) to create a mutated cycle 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 T generations or if the performance index is not improved more than ε within t s generations.

Examples
Example 1.Let the parallel chemical reactions A i → B i , A i → C i take place in two continuously stirred tank reactors with the inventory interactions (see Figure 1), where A i is the raw material of the ith reactor, B i is its desired product, and C i is its waste product (i = 1, 2).Assume x 12 (t) x 22 (t) that the ith reactor is τ i -periodically operated, and denote by x i1 (t), x i2 (t), x i3 (t) its concentrations of A i , B i , C i , respectively, by u i1 (t) its input concentration of A i , by u i2 (t) its flow rate, by v i (t) its inventory interaction transferring the waste product of the cooperating subsystem as the supplement of the raw material of the ith subsystem.Consider the following GOMC problem for the discussed system: minimize the global objective function subject for i = 1, 2 to the constraints where x i j (0) = x i j (τ i ) (j = 1, 2, 3), x i (t) .= (x i1 (t), x i2 (t), x i3 (t)) T is the state of the ith subsystem, and u i (t) .
= (u i1 (t), u i2 (t)) T is its control.The reactions obey the power law with the exponents p i j .The optimization goal is equivalent to the maximization of the summary gain of the cooperating systems for the gain coefficients c i connected with the ith useful product B i (i = 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 c i = 1, n i = 3, m i = 2, p i1 = 2, p i2 = 1 (i = 1, 2), κ 11 = 2, κ 12 = 0.5, κ 21 = 10, κ 22 = 3 we obtained the π-curves shown in Figure 2. The presented π-curves were calculated for optimal steady-state solution u 1 = u 2 = (1, 1) T , v 1 = 0.64, v 2 = 0.30, x 1 = (0.62, 0.73, 0.30) T , x 2 = (0.21, 0.45, 0.64) T which gives the performance index G = −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.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.
International Journal of Chemical Engineering   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 v 1 (t) of the first subsystem transfers the waste product of the cooperating subsystem as the catalyst for the chemical reaction where x i j (0) = x i j (τ i ) (j = 1, 2, 3).The reactions obey the power law with the exponents p i j (i = 1, 2; j = 1, 2).The results presented below were obtained for the following subsystem parameters: Just like in the previous example the optimal static control process was compared with periodic, and multiperiodic control processes.The optimal steady-state solution u 1 = u 2 = (1, 1) T , v 1 = 0.30, v 2 = 0.40, x 1 = (0.50, 0.09, 0.40) T , x 2 = (0.15, 0.94, 0.30) T gives the performance index G = −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 G = −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 τ 1 = 19.7,τ 2 = 6.5 gave the performance index G = −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.

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 finitedimensional 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.

Step 2 .
Create function w(χ k ) (k = 0, 1, . . ., K − 1), which returns values from the range 1, 2, . . ., K depending on the value of χ k .If (1/K) k∈{0,1,...,K−1} χ k ≥ γ the function w(χ k ) returns the position of element χ k in the set {χ k } K−1 k=0 sorted in nonincreasing order.Otherwise the function w(χ k ) returns the position of element χ k in the set {χ k } (i) draw from the current population β individuals and (using the formula (25)) compute for them the modified performance index G( z ν ) (ii) take from β individuals to a new ( t + 1) population one individual with the smallest value of the G( z ν ) (iii) remove from the current population β individuals.

Figure 1 :
Figure 1: Diagram of the parallel chemical reactions from Example 1.

Figure 2 :
Figure 2: The results of the π-test for Example 1.

Figure 8 :
Figure 8: The results of the π-test for Example 2.
respectively, while s i .= (s i j ) ni j=1 is the vector of the eigenvalues of the monodromy matrix Φ i (τ i , x i , u i , v i ) of the state equation (2), |s i | ∞ .= max j=1,2,...,ni |s i j |, and α i ∈ R + is the local stability level of the ith subsystem.
Figure 12: The suboptimal multiperiodic inventory interaction v 2 for Example 2.Figure 13: The suboptimal multiperiodic state x 1 for Example 2.
Figure 14: The suboptimal multiperiodic state x 2 for Example 2.