Introduction

The design of mathematical models and the development of methods for their research are necessary for an adequate description of the functioning of the new classes of technical systems. The problem of optimal parameters search for the controlled technical systems motion is an actual problem arising in the modeling of many classes of the controllable technical systems. The questions of mathematical modeling for some classes of such systems are considered in [Masina, 2006], [Masina & Druzhinina, 2016], [Shestakov et al., 2014], [Druzhinina et al., 2017], [Druzhinina et al, 2017]. The examples of such systems are aircraft control systems, space vehicles, autocontrolled cars, compact unmanned flying vehicles including multicopters. The problem of optimal control becomes especially relevant for the model of the control of several unmanned vehicles so called swarm.

In present paper we use the theory of differential inclusions and optimal control theory for description of the processes of dynamics and control of the technical systems. The practical application of the theory of differential inclusions is connected with the possibility of combining the results of this theory with the Pontryagin maximum principle, with stability theory and with methods of dynamic programming [Pontryagin, 1985], [Blagodatskikh & Filippov, 1985], [Aleksandrov et al., 2005], [Tolstonogov, 1986],

[ Borisovich et al., 2011], [Merenkov, 2000], [Shestakov, 1990], [Shestakov & Druzhinina, 2002] [Masina, 2008], [Druzhinina & Masina, 2009]. In the case when we describe the model using differential inclusions, it is convenient to write assumptions on the controlled system using the multivalued function f (x, u), where x is the trajectory of the system, u is the control function. A wide class of objects, including control systems with variable control region, as well as control systems with state constraints [Pontryagin, 1985], [Blagodatskikh & Filippov, 1985], [Aleksandrov et al., 2005] can be described using the differential inclusions.

The existence and stability of solutions of differential inclusions were considered in [Merenkov, 2000], [Tolstonogov, 1986], [Borisovich et al., 2011], [Masina, 2008] and in other works.

The synthesis and the sequential analysis of different models of controlled systems described by differential inclusions are performed in the present paper. The mathematical models of optimal control are constructed and investigated taking into account the tracking of the point moving on the plane. First, we consider a simple two-dimensional model without taking into account the variability. Secondly, a three-dimensional model is constructed and studied taking into account the variability. Thirdly, a three-dimensional model taking into account the aerodynamic forces is constructed and studied. The algorithms of the optimal parameters search for three-dimensional models are proposed. The software package that allows selecting optimal parameters of the system motion in case of changing input data is developed.

The Optimal Parameters Search of the Two-Dimensional Model

Initially, we consider the two-dimensional model of the controlled dynamical system with Cartesian coordinates x, y and starting point (0, 0). Let the object motion is carried out in two stages: in the time intervals (0, t 1 ) and (t 1 , t 2 ) for the first and second stages respectively. At the first stage the object with mass m moves in the plane xoy under the influence of the constant vector thrust (p, q) until the maximum height h. At the second step the object moves under the influence of vector thrust (b, s), the boundary conditions have the form

x(t 2 ) = l, ẋ(t 2 ) = 0, y(t 2 ) = 0, ẏ(t 2 ) = 0. Selecting of vector thrust values is made from strictly positive intervals p ∈ (p 1 , p 2 ), q ∈ (q 1 , q 2 ), b ∈ (b 1 , b 2 ), s ∈ (s 1 , s 2 ).

The considered model takes into account the impact of the gravitational force with the magnitude of the acceleration of gravity g.

The differential inclusions describing the system have the form:

{ mẍ ∈ p, mÿ ∈ q − mg, 0 < t ≤ t 1 ,(1){ mẍ ∈ −b, mÿ ∈ s − mg, t 1 < t ≤ t 2 .

(2)

The single-valued implementations of the model ( 1), ( 2) can be represented in the form of ordinary differential equations systems

{ mẍ = p, mÿ = q − mg, 0 < t ≤ t 1 ,(3){ mẍ = −b, mÿ = s − mg, t 1 < t ≤ t 2 . (4)

The optimality criterion is

∫ t1 0 (p + q)dt + ∫ t2 t1 (b + s)dt → min . (5)

From ( 3), (4) taking into account the optimality criterion (5) we obtained that the parameters p, q, b, s, t 1 , t 2 in case of single-valued realization have the form

p = 2lτ (τ − 1) 2 gm 4l(τ − 1) 2 + 2hτ 2 , q = 2hτ 2 (τ − 1) 2 gm 4l(τ − 1) 2 + 2hτ 2 + mg, b = 2lτ (τ − 1)gm 4l(τ − 1) + 2hτ 2 , s = 2hτ 2 (2τ − 1)gm 4l(τ − 1) 2 2hτ 2 + mg, τ = 1 + 1 √ 3 , t 1 = √ 4l(τ − 1) 2 + 2hτ 2 τ 2 (τ − 1) 2 g , t 2 = τ t 1 . (6)

The multivaluence that can be caused by the rarefied medium resistance factor in the described model takes into account when constructing a set of the single-valued realizations of the form (3), (4) in case when these implementations describe the motion of the object within a circle. The graph images of the possible trajectories of the object motion are presented in [Masina, 2006] when the endpoint is localized in the circle of positions of the plane xoy. The stability analysis of the nominal motion determined by the criterion of optimality is carried out [Masina, 2006].

The computational package Mathematica is used to find the formulas (6). Model (1), ( 2) is a strongly simplified model for describing the motion of the controlled object since it does not take into account the possibility of considering a mobile end point (variability condition). Besides, model ( 1), (2) does not take into account air resistance and requires generalization and clarification.

The Synthesis of the Three-Dimensional Dynamic Model

The two-dimensional model modification is considered in case when the endpoint moves in the additional plane xoz. The algorithm for constructing the motion trajectory taking into account the optimal conditions, which consists of three steps is developed.

Step 1. Object moves in the plane xoz under the influence of a constant thrust vector before reaching the point (0, h, 0). During the motion time t 1 end point moves to the point of the positions circle.

Step 2. Out of the all possible trajectories of the object motion we select the one that leads to an intermediate destination point. For the remaining time the object will reach the plane xoz taking into account the direction to the center of the circle of positions. Next we consider the motion for the time t/2 and establish a new circle of positions and select the trajectory leading to its center.

Step 3. If the radius of the positions circle is larger than the radius of the conformity, we repeat the step 2. Wherein the time of the motion, when the object passes half of the route in the center of the circle of positions is halved. As a result of several repetitions of step 2, we get a circle of positions with a radius of less than σ that allows to continue to motion in the center of the circle of positions with a guaranteed hit at the endpoint.

The two-dimensional model is modified for the implementation of the algorithm of the dynamical system motion. Let us choose the radius of the projection endpoint r and the radius conformity σ. In step 2 a sequential solution of equations system under condition of endpoint t 2 motion along axis Ox is carried out. The system describing the object motion is converted to the form

ẍ = −R 1 , x(t 1 ) = a 1 , ẋ(t 1 ) = c 1 , ÿ = −S 1 , y(t 1 ) = a 2 , ẏ(t 1 ) = c 2 , z = Z 1 , z(t 1 ) = a 3 , ż(t 1 ) = c 3 , (7)

where R 1 , S 1 , Z 1 are the acceleration components acting on the x, y, z axis, a 1 , a 2 , a 3 are the final values of previous iterations. The derivatives c 1 , c 2 , c 3 can be approximately calculated by the formula

c i = x k − x k−1 t k − t k−1 ,

where k is the number of iteration of the algorithm. Given the conditions ẋ(t 2 ) = 0, ẏ(t 2 ) = 0, ż(t 2 ) = δ, a 1 = h 2 , t 2 = τ t 1 from (7) we obtain:

R 1 = − c 1 t 1 (1 − τ ) , S 1 = − c 2 t 1 (1 − τ ) , Z 1 = − c 3 t 1 (τ − 1) + δ.

The trajectory of motion considering optimal values R 1 , S 1 , Z 1 corresponding to the single-valued implementation under conditions l = 10, h = 5, g = 9.8, m = 1, σ = 1, a 1 = 0, a 2 = 0, a 3 = 0, c 1 = 0, c 2 = 0, c 3 = 0 has the view shown in Figure 1. The set of single-valued realizations of the model describes the motion of the system taking into account the optimality and variability. Thus, the model of the dynamical system described by differential inclusions is constructed taking into account optimality and variability. The software package in the built-in language of the Octave 4.0 system is elaborated on the basis of the algorithm [Masina, 2006].

Generalization of the Three-Dimensional Model and the Optimal Parameters Search

Next, we will consider the modification of the three-dimensional model taking into account the influence of aerodynamic forces. Let the object of control has the initial coordinates (0, y 0 , 0); the initial velocity is ( x0 , ȳ0 , 0). The motion is being considered at the finite time interval T ∈ [t 0 , t ∆ ]. The endpoint moves according to a certain law, so that its position can only be determined in each moments t i , (i = 1, 2, ..., n) of motion. The center of the object mass is affected by influencing the thrust vector A as well as the force of gravity. The aerodynamic force acting on the object consists of three components -the drag, the lifting force, directed orthogonally to the motion vector, and the aerodynamic moment. The value of the axial position (α, β, γ) is also considered. Let us assume that the motion is accessible if the coordinate y is decreasing continuously and equals zero at the end point L(l x , l z , 0), as well as stable if the object axis coincides with the direction of motion. Object motion equations have the form

ẍ = 1 m (A x + (ϕ x (α, β, γ, ẋ) + ψ x (α, β, γ, ẏ, ż)) ẋ2 ), ÿ = 1 m (A y − g + (ϕ x (α, β, γ, ẏ) + ψ y (α, β, γ, ẋ, ż)) ẏ2 ), z = 1 m (A z + (ϕ x (α, β, γ, ż) + ψ z (α, β, γ, ẋ, ẏ)) ż2 ), α = 1 κ α (B α + µ α (α, β, γ, ẋ, ẏ, ż)), β = 1 κ β (B α + µ β (α, β, γ, ẋ, ẏ, ż)), γ = 1 κ γ (B γ + µ γ (α, β, γ, ẋ, ẏ, ż)),(8)

where ϕ, ψ, µ are the drag parameters, the lifting force and aerodynamic moment, respectively defined for each object form experimentally, A, B are the controlled actions in the respective measurements, κ is the axial moment of inertia. Optimality criterion is

∫ t∆ t0 (|A x | + |A y | + |A z | + |B α | + |B β | + |B γ |)dt → min .

From the algorithm in the Section 3 let us represent the time interval T as a sum of subintervals species

T = n ∑ i=0 t ∆ − t i 2 .

The final trajectory E of system can be represented as a sum of intermediate trajectories in each interval [t i , t i+1 ] ∈ T as

E = n ∑ i=0 M i (t), t i ≤ t < t i+1 ,

where M i is vector (x i , y i , z i ). The considered restrictions have the form

x i (t ∆ ) = l i , y i (t ∆ ) = 0, z i (t ∆ ) = k i , Ṁi (t ∆ ) = 0,

where l i and k i are the results of measurement of endpoint displacement.

The search for optimal values of A i and B i is the search minimization functional min(A i , B i ) = = f (P i , V i , t i , t ∆ , ϕ, ψ, µ) where P i , V i are initial coordinates and velocity vectors respectively. The algorithm of optimal parameters finding using artificial neural networks [Haykin, 2006] and logical controllers [Gostev, 2008] for model ( 8) is designed. The graph of the algorithm of the optimal parameters research includes components such as a fuzzyfication, a rule base, a defuzzyfication, a neural network, sets of input and output parameters of the model ( 8).

The analysis of the nonlinear model to assess training error criterion is performed. The training process means weighting coefficients tuning. The training criterion is the following

g(t i )     l i 0 k i   − M i (t ∆ )   2 + z(t i )     0 0 0   − Ṁi (t ∆ )   2 → min,

where M i , Ṁi calculated along the optimal values A i , B i . The functions g(t) and z(t) in each interval [t i , t i+1 ] are analogues to the membership functions, which allows us to carry out braking when approaching the endpoint. The advantage of the proposed algorithm is the opportunity to study the model in case when analytical research is difficult. Besides, the algorithm is characterized by the insignificant computational complexity. In the considered case it is possible to use Takagi-Sugeno controllers [Tanaka & Wang, 2001], [Druzhinina & Masina, 2009], which allow to carry out the "downward" method of modeling. The adjustment of the parameters with such modeling occurs in the course of a more detailed study of the dynamic system.

Conclusions

The considered models show different approaches to the search of optimal parameters of control for technical systems. The first approach is connected with construction of the dynamical model in two-dimensional space taking into account the multivalence and criteria of an optimality. The second approach is the development of the first approach and second approach is based on the construction of model in three-dimensional space taking into account the multivalence and variability. The third approach is the generalization of the first and second approaches in case when air resistance is taken into account. The algorithms for finding of the optimal parameters of control including on the basis of artificial neural networks are proposed. These algorithms can be implemented as a computer programs. The obtained results can be used for a wide range of applications, in particular, the problems of autopiloting car for the destinations taking into account topography and external conditions, the problems of autonomous control for reusable spacecraft in long-term research missions with regard to optimality of motion in gravitational complex configuration fields.