Application of the modified self-organizing migration algorithm MSOMA in optimal open-loop control problems of switching systems

. The article discusses the application of the metaheuristic algorithm of global constrained optimization for solving the problem of ﬁnding the optimal open-loop control for nonlinear switching deterministic dynamical systems. The quality of control is assessed by the value of the functional deﬁned on individual trajectories. The optimal control problem is reduced to a parametric optimization problem, which is solved using the MSOMA algorithm, which belongs to the evolutionary group. The MSOMA algorithm is a new algorithm based on the SOMA self-organizing migration algorithm. The modiﬁcation consists in identifying three leaders among the individuals forming the current population. For each member of the population, two clones are generated with the same position. Thus, in fact, three populations are generated, each of which then realizes a migration cycle (evolutionary process) relative to one of the three selected leaders. A step-by-step algorithm for piecewise-constant, piecewise-linear, quadratic spline and cubic spline methods of control laws approximation is proposed. The e ﬀ ectiveness of the proposed method is demonstrated by the example of solving optimal control problems for switching systems with two and three subsystems. The inﬂuence of the parameters of the MSOMA algo-rithm on the quality of the obtained result is investigated. Comparison of the operation of the algorithm with the known solution is carried out.


Introduction
The solution of the problem of finding the open-loop optimal control of deterministic dynamical switching systems is considered.The object model is described by a system of ordinary differential equations.The controls are constrained by a parallelepiped type.The quality of control is assessed by the value of the quality functional.
For the numerical solution of the problem of finding the optimal open-loop control of deterministic switching systems, as a rule, the necessary optimality conditions in the form of the maximum principle are applied together with methods for solving two-point boundary value problems (shooting method, residual minimization, grid method, differential sweep method, finite element method, etc.) [1].The collocation method, methods of approximation of control by piecewise-constant, piecewise-linear functions, splines, expansions in various systems of basis functions, spectral and quasi-spectral methods still remain popular [2][3][4][5][6].
In most approaches, the problem of finding the optimal control is reduced to the problem of parametric optimization, for the solution of which both classical methods of zero, first and second orders are used, as well as new metaheuristic algorithms [7,8].
In this article, to solve optimal control problems, the transition to a discrete problem is carried out, and then the solution to the original problem is constructed by interpolating the values at the grid nodes.When searching for a control in the form of a piecewise constant, piecewise linear function or splines of different orders, the calculation formulas of the Runge-Kutta method of the fourth order are used, into which expressions for the control law are substituted corresponding to the used type of approximation.Next, a single-criterion optimization algorithm is used to select the optimal values of the parameters that define the desired control.
As an optimization method, it is proposed to use a modified self-organizing migration algorithm MSOMA [9].The MSOMA algorithm is a new algorithm based on the SOMA self-organizing migration algorithm [10].A computational procedure is proposed, the effectiveness of which is demonstrated by the example of optimal control problems for a chemical process in a mixing reactor and a singular optimal control problem [7].

Statement of the problem
The behavior of the control object model is described by a system of differential equations: where x i(t) is system state vector, x i(t) = x 1,i(t) , . . ., x n,i(t) T ∈ R n ; t is time, t ∈ T = [t 0 , t f ] is the time interval of the system operation; the start t 0 and end points t f of the process are set; h is time lag; u i(t) is control vector, u i(t) = u 1,i(t) , . . ., u q,i(t) ∈ U i(t) (t) ⊆ R q ; U i(t) (t) is the set of possible control values in the mode i(t), described by the direct product of segments [a j,i(t) (t), b j,i(t) (t)], j = 1, . . ., q; f i(t) (t, x, y, u) is continuously differentiable vector function of the right parts of the subsystem equation with number i, f i(t) (t, x, y, u): The change of one subsystem to another is determined by a switching sequence σ over a period of time [t 0 , t f ]: where 0 Each element of the sequence contains the moment of switching from one subsystem to another and the number of the next (active) subsystem.Initial conditions: where x 0,i(t 0 ) is initial state of the system, ϕ i(t 0 ) (t): R 1 → R n is a piecewise continuous function of the system's prehistory.At the moments of switching, the state vector of the system generally changes according to the transition equation: where M i j is transition matrix of size (n × n), T i j ∈ R n is the transition vector from the subsystem with the number j to a subsystem with a number i.If T i j ≡ θ, where θ is zero vector, then linear jumps of the form x(t k ) = M i(t k )i(t k −0) x(t k − 0) are observed in the system, and if M i j ≡ E is a single matrix of order n, then transitions between subsystems occur while maintaining continuity of trajectories.Since the system (1) is a system with a delay, after the transition to a new subsystem, the background of the system with the number is used i(t k − 0): The set of acceptable controls U 0 form multi-link piecewise continuous functions The switching system is a corte S = {D, F}, where D(I, Γ) is a directed graph with a set of vertices I = {1, 2, . . ., M} and lots of arcs Γ, where e i, j ∈ Γ, e i, j = �i, j� is the arc coming out of the vertex.i ∈ I. and entering the vertex j ∈ I, as well as many vector functions Set of vertex heirs i this is a set of nodes reachable from a vertex i: succ(i) ∈ I, i k+1 ∈ succ(i k ), k = 0, . . ., K. The graph is given by the adjacency matrix.
The performance index to be minimized is given by where � is penalty function for switching from one subsystem with a number i(t k − 0) to another subsystem with a number i(t k ).
In order to calculate the value of the functional I it is necessary to find the trajectory of the system x(t) = � i K i(t)=i 0 x i(t) (t), formed by the trajectories of subsystems (links) determined by the switching sequence (2) corresponding to the permissible control u(t), from the equation of the system (1) taking into account the initial condition (3) and the transition equation (4).
It is required to find such a switching sequence σ * = � (t 0 , i 0 ), (t 1 , i 1 ), . . ., (t from the set of permissible controls on which the minimum value of the functional I is achieved.

Solution search strategy
The basis of the strategy is the decomposition of control in the interval [t k , t k+1 ) according to the basis system of functions [11] using the saturation function, taking into account the specified control constraints: where sat ), j = 1, . . ., q, N is the number of unknown decomposition coefficients set by the user.
At the same time, the structure of the block matrix is a column of optimized control parameters The optimization problem is to choose the vector coefficients by minimization of the cost functional value (5).

Algorithm
Step 1. Set: K max is a maximum number of switches from one subsystem to another; N max is a maximum number of intervals of the sign of the control of each subsystem; N K is a maximum number of attempts to improve the number of switches; N sw is a maximum number of attempts to improve the switching moments from one subsystem to another with a fixed number of switches; N mode is a maximum number of attempts to search for sequences of changing mode numbers with a fixed number of switches and switching moments; N contr is a maximum number of attempts to search for subsystem control values.
Step 2. Start the iteration counter to improve the number of switches.
Step 3. (The first level of the search procedure).Generate a random number of switches in a switching sequence K ∈ {0, 1, . . ., K max } and fix it.One can generate a number K, using a uniform distribution law on a segment [0, K max ], rounding with a disadvantage, removing the value already used in the search from the set of acceptable values.One can also start with K = 0, by increasing the number of switches by one each time Check the iteration counter to improve the number of switches.If the number N K is reached, then complete the procedure and proceed to step 9.If not, then enter the resulting solution into the Pool set and continue.
Start the counter of iterations of the search for switching moments.
Step Start the counter of attempts to search for a sequence of subsystem shift numbers.
Step 5. (The third level of the search procedure.)Randomly generate a sequence based on a given adjacency matrix σ = � (t 0 , i 0 ), (t 1 , i 1 ), . . ., (t K , i K ) � , containing information about the change of one subsystem to another.From the set of heirs of each of the vertices of the graph, one of the possible heirs is selected: i 0 ∈ I, i k+1 ∈ succ(i k ), i k+1 ∈ I, k = 0, 1, . . ., K − 1. Fix the resulting sequence.Check the iteration counter of searching for a sequence of subsystem shift numbers.If the number N mode is reached, then go to step 4.
Start the counter for the number of attempts to search for subsystem control values.
Step 6. (The fourth level of the search procedure).
Set the control on the interval [t k , t k+1 ) in the form of a decomposition of the desired control law according to a basic system of functions using a saturation function that takes into account the specified control constraints (6).
Thus, the control of the switching system can be represented by a block matrix-column (7) Step 7. Using the control u(t) = � i K i(t)=i 0 u i(t) (t), set by the matrix A, find the trajectory of the switching system x(t) = � i K i(t)=i 0 x i(t) (t) and calculate the value of the control quality functional ( 5): To calculate the functional, it is required to integrate the system (1) together with the equation for the auxiliary variable: ẋi ),n+1 (t k ) using the switching sequence (2), the initial conditions (3) and the transition equation ( 4) using the MSOMA [12] method.
Step 8. Go to the next matrix-column to find the minimum functional I by one of the metaheuristic methods, generating a new piecewise constant control for each of the active subsystems (step 6).Check the counter of the number of attempts to search for subsystem control values.If the number N contr is reached, then go to step 5.
Step 9. Choose the best solution from the Pool set.

Numerical examples Example 1
Consider a switching system consisting of three nonlinear subsystems in table 1.
Let t 0 = 0, t f = 3.The system switches at a moment in time t = t 1 from subsystem 1 to subsystem 2 and at a time t = t 2 from subsystem 2 to subsystem 3 (0 ≤ t 1 ≤ t 2 ≤ 3).Initial conditions: x 1 (0) = 2, x 2 (0) = 3.Initial switching time values: It is required to find optimal switching moments t 1 , t 2 and optimal control u(t), providing a minimum of the cost functional: Table 2 shows the results of the influence analysis of the MSOMA parameters for various types of control law approximations.The known solution was obtained in [13]: t opt 1 = 0.2262, t opt 2 = 1.0176,I opt = 5.4399.So, the best solution obtained by the presented algorithm is better by about 0.005.

Example 2
Consider a switching system consisting of two linear subsystems in table 3.
Let t 0 = 0, t f = 2.The system switches once at a time t 1 from subsystem 1 to subsystem 2 (0 ≤ t 1 ≤ 2).In this case, switching occurs with a discontinuous jump: It is required to find the optimal switching moment t 1 and optimal control u(t), providing a minimum of the cost functional: Table 4 shows the results of the influence analysis of the MSOMA parameters for various types of control law approximations.In this work, a strategy, a step-by-step algorithm and corresponding software for the approximate solution of the problem of optimal open-loop control of switching systems have been developed.The presented algorithm and program are tested on the examples of solving optimal control problems for switching systems with two and three subsystems.The influence of the parameters of the MSOMA algorithm on the quality of the obtained result is investigated.Comparison of the operation of the algorithm with a known solution, as well as with various methods of control laws approximation is given.Recommendations on the choice of algorithm parameters are given.

Figure 2 .
Figure 2. Control law using piecewise linear approximation (a) and phase trajectory (b) Check the iteration counter of the search for switching moments.If the number N sw is reached, then go to step 3.
4. (Second level of the search procedure.)The time interval of the system operation [t 0 , t f ] should be divided into segments [t 0 , t 1 ], [t 1 , t 2 ], . . ., [t K , t f ], where {t k } K k=1 is a sequence of switching moments generated randomly with a condition check t 0 ≤ t 1 ≤ . . .≤ t K ≤ t f ; fix

Table 2 .
Influence of MSOMA parameters for different types of approximation

Table 4 .
Influence of MSOMA parameters for different types of approximation