# Solving Boundary Value Problem for a Nonlinear Stationary Controllable System with Synthesizing Control

Authors: Alexander N. Kvitko; Oksana S. Firyulina; Alexey S. Eremin
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

An algorithm for constructing a control function that transfers a wide class of stationary nonlinear systems of ordinary differential equations from an initial state to a final state under certain control restrictions is proposed. The algorithm is designed to be convenient for numerical implementation. A constructive criterion of the desired transfer possibility is presented. The problem of an interorbital flight is considered as a test example and it is simulated numerically with the presented method.


---

## Body

## 1. Introduction and Problem Formulation

Transferring an object from a known initial state to a desired point in the phase space is an important problem of the motion control theory. It is quite difficult and has not been solved in the general case. Two classic approaches to solve it are program trajectories and targeting trajectories. In the former case, the object moves along the fixed trajectory which cannot be changed (i.e., the problem of its stabilization arises). This leads to strong restrictions on the control system possibilities. In the latter case, the motion is controlled constantly (with some time step) based on the current object state in order to perform the desired transfer. This demands constant target point tracking with some technical means, which reduces the control system independence.Nowadays, the progress in computers, measurement tools, and other instruments’ engineering makes it possible to start development and production of autonomous intellectual control systems, which themselves can determine the transfer trajectory with onboard computational and navigational systems using only the initial (or current) state of the object. This greatly extends the capabilities of control systems, since there is no need in constant target point tracking, which often cannot be practically realized.Intellectual control systems for many technical objects (aerospace vehicles, gyroscopic mechanisms, robotic manipulators, etc.) are modeled mathematically with nonlinear stationary controllable systems of ordinary differential equations (ODEs). The necessary control action is found as control functions, which provide the notion that phase coordinates satisfy the given boundary conditions. Such problems of the mathematical control theory are called boundary value problems (BVPs) for ODEs. They form an important class of problems and are quite difficult, especially when additional details are introduced into the model, like discrete time consideration, control signal and measurement delays, permanent perturbations, or integral-differential components in the object behavior.For the first time, the complete solution of the boundary value problem for linear controllable systems in a class of square-integrable control functions was proposed by Kalman et al. [1]. During the following decades, a lot of papers and monographies studying BVPs for linear and nonlinear controlled systems of ODEs and integral-differential equations of a special kind have appeared [1–33]. It should be also mentioned that the control methods based on the classic approaches continue to appear, especially in the areas where models become even more complex, as fuzzy control systems (see, e.g., [34–37]).Today, there are three main directions of BVP study. The first deals with finding the necessary and sufficient conditions imposed to the right part of the controlled system which guarantee the object transfer to the desired point in the phase space [1–5, 7–10, 13–16, 18, 19, 21–32].The second problem is to determine the final states set, where the object can be transferred from the initial state [3, 5–8, 11, 17, 20–22, 33].The third direction concerns development of exact or approximate methods for constructing or synthesizing program control functions and the corresponding trajectories which connect the given points in the phase space [1, 3–6, 11, 18, 20–22, 24, 26].BVPs’ solution is studied sufficiently well for linear controllable systems. For quasi-linear (semilinear) and nonlinear systems of the special forms, mainly iterative methods and their modifications are used [4, 5]. The weaknesses of such methods are instability to computational errors and outer perturbations and high computational time, since each stage includes the solution of a nonlinear system of functional equations. Moreover, it is quite difficult to find good initial approximations to the control functions and time functions of phase coordinates to provide fast convergence. Initial approximations are often chosen intuitively. Due to the mentioned problems, it is impossible to use successive approximations to solve BVPs in essentially nonlinear controllable systems. Theory and methods of solving BVPs for nonlinear systems of general form have yet to be developed.The present paper is based on the results presented in [2, 24]. In [2], the condition of Kalman type, which guarantees the local controllability for a wide class of nonlinear stationary systems of ODEs, is considered. In [24], the problem of local controllability of a nonlinear stationary system is considered. The notions of auxiliary controllable systems linearized along the trajectory xt and control ut and auxiliary controllable system linearized about the original nonlinear system equilibrium point xε,uε are introduced. It is proven that the local controllability of the original nonlinear system follows from the full controllability of the auxiliary systems.In the current paper, a nonlinear stationary system local controllability condition is found, which is analogous to the condition presented in [2, 24]. The main difference from the previous results is that we develop an algorithm for constructing the synthesizing control function, which provides the transfer from the origin of coordinates to a given final state for a wide class of nonlinear stationary systems of ODEs under restrictions on the control action magnitude. The algorithm is quite easy to implement and it has stable input and computational errors. Moreover, the estimations of the reachable final states set are presented. The algorithm realization is simple: the original problem is reduced to the stabilization of a nonlinear system of the special form under permanent perturbations, and then an initial value problem (IVP) for an auxiliary system of ODEs is solved. The most complex and time-consuming part, which is the auxiliary system construction and its stabilization, can be made with analytical methods and realized with computer algebra means. Additionally, the required control law found from the auxiliary system stabilization condition provides the stability of computational procedures to input and computational errors. The effectiveness of the presented algorithm is demonstrated with numerical solution of a practical problem.In the paper, we consider a controllable system of ordinary differential equations (ODEs):(1)x˙=fx,u,where x=x1,…,xnT is a vector containing n phase coordinates and control u is a vector of the same or lesser dimension: u=u1,…,urT and r≤n. We choose the time scale so that the considered time t∈[0,1]. The right-hand side is(2)f∈C4nRn×Rr;Rn,f=f1,…,fnT,(3)f0,0=0.Notice that(4)rank⁡S=n,where(5)A=∂f∂x0,0,B=∂f∂u0,0,S=B,AB,A2B,…,An-1B.The possible control magnitude is bounded:(6)u<N.Problem 1.
Find a pair of functionsxt∈C0,1 and ut∈C0,1 that satisfy (1) and the boundary conditions(7)x0=0,x1=x-,x-=x-1,…,x-nT.We say that the pair xt and ut is a solution of problems (1) and (7).Theorem 2.
If conditions (2), (3), and (4) are satisfied for the right-hand side of (1), then ∃ε>0 such that ∀x-∈Rn:x-<ε there exists a solution of problems (1) and (7), which can be found after solving, first, a problem of stabilizing a linear nonstationary system with exponential coefficients and, second, an initial value problem for an auxiliary ODE system.The main idea of the proof is to use successive changes of independent and dependent variables to reduce the process of solving the original system to the problem of stabilizing a nonlinear auxiliary system of ODEs of the special form under permanent perturbations. To solve the latter, we find a synthesizing control, which provides exponential decrease of the linear auxiliary system fundamental matrix. At the final stage, we return to the original variables.


## 2. Auxiliary System Construction

We find the functionxt being a part of the solution of (1) and (7) in the form(8)xit=ait+x-i,i=1,…,n.In the new variables, system (1) and the boundary conditions are written as (9)a=fx-+a,u,(10)a0=-x-,a1=0.We call the pair of functionsat and ut, which satisfy system (9) and conditions (10), a solution of (9) and (10). Now consider a problem of finding at∈C0,1 and ut∈C0,1, which satisfy (9), such that(11)a0=-x-,at⟶0ast⟶1.Remark 3.
The limit witht→1 of the solution of (9) and (11) is the solution of (9) and (10).We switch from the independent variablet in system (9) to τ according to(12)t=1-e-ατ,τ∈0,+∞,where α>0 is a certain constant value to be determined. Then, in terms of τ, (9) and (11) take the form(13)dcdτ=αe-ατfx-+c,d,τ∈0,+∞,(14)c0=-x-,cτ⟶0asτ⟶∞,(15)cτ=atτ,dτ=utτ,c=c1,…,cnT,d=d1,…,drT.We call the pair of functions cτ and dτ, which satisfy system (13) with conditions (14), a solution of problem (13) and (14). With the solution of (13) and (14), one can return to the solution of (9) and (11) with (12) and (15).Let us denote(16)c~=x-+θic,d~=θid,θ∈0,1,i=1,…,n,k=∑i=1nki,m=∑i=1rmi,k!=k1!·…·kn!,m!=m1!·…·mr!.Using property (2) and Taylor series expansion of the right-hand side of (1) about x-,0, we can rewrite system (13) as(17)dcidτ=αe-ατfix-,0+αe-ατ∑j=1n∂fi∂xjx-,0cj+αe-ατ∑j=1r∂fi∂ujx-,0dj+12αe-ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0cjck+2∑j=1n∑k=1r∂2fi∂xj∂ukx-,0cjdk+∑j=1r∑k=1r∂2fi∂uj∂ukx-,0djdk+αe-ατ∑k+m=4n-21k!m!∂k+mfi∂x1k1⋯∂xnkn∂u1m1⋯∂urmrx-,0c1k1×⋯×cnknd1m1×⋯×drmr+αe-ατ∑k+m=4n-11k!m!∂k+mfi∂x1k1⋯∂xnkn∂u1m1⋯∂urmrc~,d~c1k1×⋯×cnknd1m1×⋯×drmr,i=1,…,n.Let us bound the range ofcτ with(18)cτ<C1,τ∈0,∞.We will now shift the functions ciτ, i=1,…,n, a number of times. Our aim is to get an equivalent system where all the terms in the right-hand side, which do not contain powers of c or d in explicit form, would be of the order Oe-4nατx- as τ→∞ and x-→0 in domain (6) and (18).At the first stage, we changeciτ to ci1τ by the following rule:(19)ciτ=ci1-e-ατfix-,0,i=1,…,n.Let Dk+mfi≡∂k+mfi/∂x1k1⋯∂xnkn∂u1m1⋯∂urmr for all i=1,…,n. After substitution of (19) into the left- and right-hand sides of (17), we obtain the equivalent system(20)dci1dτ=-αe-2ατ∑j=1n∂fi∂xjx-,0fjx-,0+12αe-3ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0fjx-,0fkx-,0+αe-ατ∑j=1n∂fi∂xjx-,0cj1+e-2ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0fkx-,0cj1+αe-ατ∑j=1r∂fi∂ujx-,0dk+e-2ατ∑j=1n∑k=1r∂2fi∂xj∂ukx-,0fjx-,0dk+12αe-ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0cj1ck1+αe-ατ∑j=1n∑k=1r∂2fi∂xj∂ukx-,0dkcj1+12αe-ατ∑j=1r∑k=1r∂2fi∂uj∂ukx-,0djdk+⋯+αe-ατ∑k+m=4n-21k!m!Dk+mfix-,0c11-e-ατf1x-,0k1×⋯×cn1-e-ατfnx-,0knd1m1×⋯×drmr+αe-ατ∑k+m=4n-11k!m!Dk+mfic~,d~c11-e-ατf1x-,0k1×⋯×cn1-e-ατfnx-,0knd1m1×⋯×drmr,i=1,…,n.It follows from (14) and (19) that(21)ci10=-x-i+fix-,0,i=1,…,n.It is easy to see that, in the right-hand side of (20), the terms, which do not contain the powers of the components of c or d in explicit form, are bounded with Oe-2ατx- as τ→∞ and x-→0 in (6) and (18).At the second stage, we make a change of variables(22)ci1τ=ci2τ+e-2ατϕi2x-,ϕi2x-=12∑j=1n∂fi∂xjx-,0fjx-,0,ϕi20=0,i=1,…,n.With respect to these new variables, the original system (20) and the initial conditions (21) are written as(23)dci2dτ=α12e-3ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0fjx-,0fkx-,0+e-3ατ∑j=1n∂fi∂xjx-,0ϕj2-e-4ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0fkx-,0ϕj2+12e-5ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0ϕj2ϕk2+αe-ατ∑j=1n∂fi∂xjx-,0cj2-e-2ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0fkx-,0cj2+e-3ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0ϕk2cj2+αe-ατ∑k=1r∂fi∂ujkx-,0dk-e-2ατ∑j=1n∑k=1r∂2fi∂xj∂ukx-,0fjx-,0dk+e-3ατ∑j=1n∑k=1r∂2fi∂xj∂ukx-,0ϕj2dk+12αe-ατ∑j=1n∑k=1n∂2fi∂xj∂xkx-,0cj2ck2+∑j=1n∑k=1r∂2fi∂xj∂ukx-,0dkcj2+∑j=1r∑k=1r∂2fi∂uj∂ux-,0djdk+αe-ατ∑k+m=4n-21k!m!Dk+mfix-,0c12+e-2ατϕ12-e-ατf1x-,0k1×⋯×cn2+e-2ατϕn2-e-ατfnx-,0knd1m1×⋯×drmr+αe-ατ∑k+m=4n-11k!m!Dk+mfic~,d~c12+e-2ατϕ12-e-ατf1x-,0k1×⋯×cn2+e-2ατϕn2-e-ατfnx-,0knd1m1×⋯×drmr,i=1,…,n.(24)c120=-x-i+fix-,0-ϕi2x-,i=1,…,n.Compared to the previous variables shift, the right-hand side terms of (23), which do not contain the powers of the components of c or d in explicit form, are now of the order Oe-3ατx- as τ→∞ and x-→0 in domain (6) and (18).By induction, at thekth stage with (19)–(24), we have the necessary shift of the form(25)cik-1τ=cik+e-kατϕikx-,ϕik0=0,i=1,…,n.We apply (25) 4n-1 times and collect the terms, which are linear with respect to the components of c(4n-1) and include the coefficients e-iατ, i=1,…,n, and also the terms, which are linear with respect to the components of d and include the coefficients e-iατ, i=1,…,2n. Now, we have the system, which according to (20)–(25) can be written in vector form as follows:(26)dc4n-1dτ=Pc4n-1+Qd+R1c4n-1,d,x-,τ+R2c4n-1,d,x-,τ+R3c4n-1,d,τ+R4x-,c4n-1,d,τ,R1=R11,…,R1nT,R2=R21,…,R2nT,R3=R31,…,R3nT,R4=R41,…,R4nT.The functionsR1i includes all the terms which are linear with respect to the components of c(4n-1) with coefficients e-iατ, i≥n+1, and also the terms of the last sum of the right-hand side for which m=0 and k=1. The function R2i includes all the terms which are linear with respect to the components of d with coefficients e-iατ, i≥2n+1, and also the terms of the last sum of the right-hand side for which m=1 and k=0. In R3i, all the terms, which are nonlinear with respect to the components of c(4n-1) or d, are contained. Finally, the function R4i includes all the terms, which do not have powers of c(4n-1) and d components. The functions P and Q have the form(27)Px-=αe-ατP1x-+e-ατP2x-+⋯+e-n-1ατPn-1x-,P1x-=∂f∂xx-,0,P10=A,Qx-=αe-ατQ1x-+e-ατQ2x-+⋯+e-2n-1ατQ2n-1x-,Q1x-=∂f∂ux-,0,Q10=B.(28)c4n-10=-x-+fx--ϕ2x--⋯-ϕn-1x-,ϕi=ϕ1i,…,ϕniT,i=1,…,4n-1,ϕi0=0.


## 3. Right-Hand Side Terms Evaluation

It follows from the construction of (26) that, in domain (6) and (18), the following estimations are true:(29)Pix-⟶0,Qjx-⟶0asx-⟶0,i=2,…,n-1,j=2,…,2n-1;(30)R1c4n-1,d,x-,τ≤e-n+1ατL1x-c4n-1,R2c4n-1,d,x-,τ≤e-2n+1ατL2x-d,L1>0,L2>0;(31)R3c4n-1,d,τ≤e-ατL3c4n-12+d2,L3>0.Moreover, conditions (2) and (3) lead to (32)R4c4n-1,d,x-,τ≤L4x-e-4nατ,L4>0.Estimation (32) follows from the representation(33)fx-,0=∂f∂xθx-,0x-,θ=θ1,…,θnT,θi∈0,1,i=1,…,n.Remark 4.
Let us denote theith column of Q1 as q1i. Construct the matrix(34)S1=q11,P1q11,…,P1k1-1q11,q12,P1q12,…,P1k2-1q12,…,q1r,…,P1kr-1q1r.Here, for each j=1,…,rkj is the maximal number of columns of the form q1j,…,P1kj-1q1j such that all the vectors q11, P1q11,…,P1k1-1q11, q12, P1q12,…,P1k2-1q12, q1r,…,P1kr-1q1r are linearly independent.It follows from (4) and (27) that there exists ε1>0 so that rank⁡S1=n for any x-∈Rn satisfying x-<ε1.We now consider the system(35)dc4n-1dτ=Pc4n-1+Qd.


## 4. Auxiliary Lemma

Lemma 5.
Let conditions (2) and (4) be satisfied for system (1). Then, ∃ε->0, ε-<ε1 such that ∀x-∈Rn:x-<ε- there exists the control dτ of the form(36)dτ=Mτc4n-1that provides an exponential decrease of the fundamental matrix in (35).Proof.
We use Krasovskii’s linear nonstationary systems stabilization method [3] in the proof. Let L1j, j=1,…,r, be the jth column of Q. Construct the matrix(37)S2=L11,L21,…,Lk11,L12,…,Lk22,…,L1r,…,Lkrr,Lij=PLi-1j-dLi-1jdτ,j=1,…,r,i=2,…,kj.Here, for each jkj is the maximal number of columns L1j,…,Lkjj such that the vectors L11, L21,…,Lk11, L12,…,Lk22, L1r,…,Lkrr are linearly independent. We show that (38)rank⁡S2=n.
LetL-1j, j=1,…,r, be the jth column of the matrix αe-ατQ1. Consider the matrix(39)S3=L-11,L-21,…,L-k11,L-12,…,L-k22,…,L-1r,…,L-krr,L-ij=αe-ατP1L-i-1j-dL-i-1jdτ,j=1,…,r,i=2,…,kj,where kj are determined the same way as for S2. Conditions (2), (27), and (29) provide the notion that S2→S3 as x-→0. This leads to the existence of ε2>0: ε2<ε1 such that ∀x-∈Rn: x-<ε2(40)rank⁡S2=rank⁡S3.Now, taking into account Remark 4 and equality (40), we can prove by contradiction that ∀x-∈Rn:x-<ε2(41)rank⁡S3=n.
From (40) and (41), it follows that condition (38) is satisfied in the domain x-<ε2. Moreover, due to the structure of S3, the estimation(42)S2-1=Oenατ,τ⟶∞,is true in that domain.
With the use of (38), we change the variable c(4n-1) according to the expression(43)c4n-1=S2τy.As a result, we get the system of ODEs:(44)dydτ=S2-1PS2-dS2dτy+S2-1Qd.
According to [38] for the first term of the right-hand side of (44) we have (45)S2-1PS2-dS2dτ=e-2,…,e-k1,φ-k1τ,…,e-k1+⋯+kr-1+2,…,e-k1+⋯+kr,φ-krτ,where e-i is a zero vector of length n with the only unit at the ith place,(46)φ-kj=-φk11,…,-φk1k1,…,-φkj1,…,-φkjkj,0,…,0n×1T,and φkji is the coefficient of Lkj+1j expansion into the sum of the vectors L11, L21,…, Lk11, L12,…,Lk22, L1r,…,Lkrr (notice that ∑j=1rkj=n); that is,(47)Lkj+1j=-∑i=1k1φk1iτLi1-⋯-∑i=1kjφkjiτLij.The second term is(48)S2-1Q=e-1,…,e-ki+1,…,e-γ+1,γ=∑i=1r-1ki.
Consider the stabilization of the system(49)dykidτ=e-2ki,…,e-kiki,φ--kiyki+e-1kidi,yki=yki1,…,ykikiki×1T,i=1,…,r,where e-iki is a zero vector of length ki with the only unit at the ith place and φ--ki=-φki1,…,-φkikiki×1T.
Letykiki=ψ. The phase variables of (49) are connected to ψτ and its derivatives with the equalities(50)ykiki=ψ,ykiki-1=ψ1+φkikiψ,ykiki-2=ψ2+φkikiψ1+dφkikidτ+φkiki-1ψ,⋮yki1=ψki-1+rki-2τψki-2+⋯+r1τψ1+r0τψ.The differentiation of the last equality in (50) together with (49) leads to(51)ψki+εki-1τψki-1+⋯+ε0τψ=di,i=1,…,r.The functions rki-2τ,…,r0τ, εki-1τ,…,ε0τ in (50) and (51) are linear combinations of the functions φkjiτ, i=1,…,kj, j=1,…,r, and their derivatives.Remark 6.
From the structure of the matricesP and Q, it follows (see (27) and representation (47)) that the functions φkikiτ,…,φki1τ and their derivatives, as well as rki-2τ,…,r0τ, εki-1τ,…,ε0τ, are bounded. Other elements of the columns φ-kj, j=1,…,r, obey the estimation Oen-1ατ, τ→∞.Let(52)di=∑j=1kiεki-jτ-γki-jψki-j,i=1,…,r,where γki-j, j=1,…,ki, are chosen to make the roots λki1,…,λkiki of the equation(53)λki+γki-1λki-1+⋯+γ0=0,i=1,…,r,satisfy the conditions (54)λkii≠λkijifi≠j,λkii<-2n+1α-1forj=1,…ki,i=1,…r.Due to (50) and (54) and Remark 6, the control di=δkiTki-1yki, i=1,…,r, provides the exponential decrease of system (44) solutions. Switching back to the original variables in (52) and using (43), we have(55)di=δkiTki-1S2ki-1c4n-1,i=1,…,r,where δki=εki-1τ-γki-1,…,ε0τ-γ0; Tki is the inequality (50) matrix, which means that yki=Tkiψ-, ψ-=ψki-1,…,ψT; S2ki-1 is the matrix which consists of the corresponding ki rows of S2-1. The desired control can be written in the form of (36), where(56)Mτ=δkTk-1S2k-1≡δk1Tk1-1S2k1-1,…,δkrTkr-1S2kr-1T.Denote the fundamental matrix of system (51) closed with control (52) via Ψτ (Ψ(0)=I, an identity matrix). It is obvious that the elements of Ψτ are exponential functions of negative argument or their derivatives.Consider system (35) with control (55):(57)dc4n-1dτ=Dτc4n-1,Dτ=Pτ+QτMτ.Introduce a block-diagonal matrix T(τ). Its diagonal blocks are matrices Tki, i=1,…,r. Then, corresponding to (43) and (50), the fundamental matrix Φτ, Φ-1(0)=I, of system (57) has the form (58)Φτ=S2τTτΨτT-10S2-10.The estimation(59)Φτ≤Ke-ατe-λτ,λ>0,K>0,follows from (58), the structure of matrices S2τ and Ψ(τ), and Remark 6.Letε-=ε2. Then, the correctness of the lemma becomes clear from (59).Besides, on base of (42), (55), and (58) and Remark 6, we have(60)ΦτΦ-1t≤K1e-λτ-ten-1αt,t≤τ,τ∈0,∞,K>0;Mτ=Oenατ,τ⟶∞.


## 5. Main Theorem Proof

System (26) with control (55) has the form(61)dc4n-1dτ=Dτc4n-1+R1c4n-1,Mτc4n-1,x-,τ+R2c4n-1,Mτc4n-1,x-,τ+R3c4n-1,Mτc4n-1,τ+R4x-,c4n-1,Mτc4n-1,τ.Change of the variables(62)c4n-1=ze-3nατ,c4n-10=z0leads to the following representation:(63)dzdτ=Cτz+e3nατR1e-3nατz,Mτe-3nατz,x-,τ+e3nατR2e-3nατz,Mτe-3nατz,x-,τ+e3nατR3e-3nατz,Mτe-3nατz,τ+e3nατR4x-,e-3nατz,Mτe-3nατz,τ,Cτ=Dτ+3nαE.Let us show that all solutions of (63) with initial values (62) that start in a sufficiently small neighborhood of zero decrease exponentially. Let Φ1τ, Φ1-10=I, be the fundamental matrix of the system dz/dτ=Cτz. Then, on the basis of (59), (60), and (62), we have(64)Φ1τ≤Ke-βτ,Φ1τΦ1-1t≤K1e-βτ-ten-1αt,β=λ-3nα.Let us choose α to provide(65)β>0.The solution of (63) with initial values (28) can be presented as(66)zτ=Φ1τΦ1-1τ1zτ1+∫τ1τΦ1τΦ1-1te3nαt×R1e-3nαtz,Mte-3nαtz,x-,t+R2e-3nαtz,Mte-3nαtz,x-,t+R3e-3nαtz,Mte-3nαtz,t+R4x-,e-3nαtz,Mte-3nαtz,tdt,forτ∈τ1,∞,(67)zτ=Φ1τc4n-10+∫0τΦ1τΦ1-1te3nαt×R1e-3nαtz,Mte-3nαtz,x-,t+R2e-3nαtz,Mte-3nαtz,x-,t+R3e-3nαtz,Mte-3nαtz,t+R4x-,e-3nαtz,Mte-3nαtz,tdt,forτ∈0,τ1.Now, from (66) and (67), using (30), (31), (32), (36), (62), and (64) after changing c to z, the following estimations in domain (6) and (18) are true:(68)zτ≤Ke-βτΦ1-1τ1zτ1+∫τ1τe-βτ-tK1L-e-αtzt+L4x-e-αtdt,forτ∈τ1,∞;zτ≤Ke-βτc4n-10+∫0τe-βτ-tK1L-e-αtzt+L4x-e-αtdt,forτ∈0,τ1,L->0.The constant L- depends on (6) and (18).We apply the well-known result from [39] to inequalities (68) and obtain(69)zτ≤Ke-μτΦ-1τ1zτ1+K1∫τ1τe-μτ-tL4x-e-αtdt,forτ∈τ1,∞,μ=β-K1L-e-ατ1;zτ≤Ke-μ1τc4n-10+K1∫0τe-μ1τ-tL4x-e-αtdt,forτ∈0,τ1,μ1=β-K1L-.Using (65), we fix τ1>0 so that the inequality μ>0 is satisfied.Now, let us bound the choice ofα>0 with the condition α<μ. Then, after the integration in the right-hand sides of (69), we have(70)zτ≤Ke-μτΦ1-1τ1zτ1+K2e-ατL4x-,forτ∈τ1,∞,zτ≤K3c4n-10+K4L4x-,forτ∈0,τ1.Note that all Ki>0.On the basis of (28) and (32), two last estimations can be written as a single inequality in the domain x-<ε-:(71)zτ≤K5e-ατx-,τ∈0,∞,K5>0.The constant K5 depends on domain (6) and (18).With the use of (36), (60), (62), and (71), we estimate dτ:(72)dτ≤Mτe-3nατK5e-ατx-≤K6e-2n-1ατx-,K6>0,forτ∈0,∞.Now, we can set the valueε from the theorem formulation, to be minε-,C1/K5,N/K6. This shows that the solution of system (63) with initial values (62) and (28) does not leave the area z<C1 for x-<ε and decreases exponentially. Besides, by (72), the corresponding function dτ obeys restriction (6). Moreover, if we substitute the solution to (62) and (36) and return to the variables cτ according to (19) and (22)–(25), we get the solution of (13) and (14).Then, we write everything in the original variables with (15) and (12) and finally (8). Now, according to Remark 3, passing to the limit as t→1 gives the solution of the original problem (1) and (7). This proves the theorem.The algorithm of solving problem (1) and (7) consists of the following steps:(1)
Construction of the auxiliary system (26); performed with analytical methods.(2)
Stabilization of system (35); the synthesizing control (55) is found; performed with analytical methods.(3)
Solution of the IVP (26) and (28) with control (55); its solution is substituted into (55) and thus the functions c4n-1τ and dτ are found; performed with numerical integration methods.(4)
Return to the original variables according to (22)–(25), (19), (15), (12), and (8); performed with analytical methods.Remark 7.
Consider the system(73)x˙=fx,u+F,where x=x1,…,xnT, x∈Rn; u=u1,…,urT, u∈Rr; t∈0,1, r≤n, F=F1,…,FnT, F∈Rn. F is a constant perturbation. From the theorem follows Corollary 8.Corollary 8.
If conditions (2)–(4) are satisfied for system (73), then there exists ε>0 such that, for any x-∈Rn and F∈Rn if x-<ε and F<ε, a solution of (73) and (7) exists an can be found after stabilizing a linear nonstationary system with exponential coefficients and solving the initial value problem for an auxiliary system of ODEs.


## 6. Example of the Algorithm Application

To demonstrate the effectiveness of the proposed method, we consider a problem of transferring a material point (a satellite) moving in a central gravitational field to a desired circular orbit with jet power.According to [3], the systems in deviations from the prescribed circular motion and condition (7) take the form(74)x˙1=x2,x˙2=-νr0+x12+χ0+x32r0+x13+aru,x˙3=x1+r0aψu,(75)xi0=0,xi1=x-i,i=1,2,3,u<N.Phase coordinates x1=r-r0, x2=r˙, and x3=χ-χ0 show the deviation from the original circular orbit with radius r0; r˙ is the radial velocity; χ is the generalized momentum and χ0=νr01/2; ar and aψ are the relative velocity vector projections onto the radial and tangential directions, respectively (they are constant); m is the mass, m˙ is its change rate and the control u=m˙/m; ν=ν0M, where ν0 is the universal gravitational constant; M is the mass of Earth. The vector of phase coordinates x=x1,x2,x3T. The control u is scalar.System (13) and conditions (14) for the problem have the form(76)dc1dτ=e-ατc2,dc2dτ=e-ατgc1,c3+e-ατard,dc3dτ=e-ατc1+x-1+r0aψd,with(77)gc1,c3=-νr0+c1+x-12+χ0+c3+x-32r0+c1+x-13,(78)ci0=-x-i,ciτ⟶0asτ⟶∞,i=1,2,3.In the following, we assume x-=x-1,0,0T. It is necessary to make seven transforms of type (25) to solve problem (76) and (78):(79)c2=z2-e-ατgx-1,c1=z1+12e-2ατgx-1,z2=w2-16e-3ατ∂g∂c1x-1gx-1,z1=w1+124e-4ατ∂g∂c1x-1gx-1,w2=υ2-e-5ατg-,w1=υ1+e-6ατg--,υ2=u2-e-7ατg---,g-=1120∂g∂c1x-12gx-1-140∂2g∂c12x-1g2x-1,g--=1720∂g∂c1x-12gx-1+1240∂2g∂c12x-1g2x-1,g---=171720∂g∂c1x-13gx-1+148∂3g∂c13x-1g3x-1+140∂g∂c1x-1∂2g∂c12x-g2x-1.Let ar=1. The matrices P and Q and the system analogous to (35) look like(80)dc-dτ=Pc-+Qd,c-=υ1,u2,c3T,P=αe-ατ010a210a23000,Q=0αe-αταe-ατb+α2e-3ατaψ,b=r0+x-1aψ,a21=∂g∂c1x-1,a23=∂g∂c3x-1.To stabilize system (80), we construct the matrix S=L1,L2,L3, with L1=Q, L2=PL1-(d/dτ)L1, L3=PL2-(d/dτ)L2; that is,(81)L1=0αe-αταe-ατb+e-2ατhL2=α2e-2ατα2e-ατe-ατa23b+e-3ατa23h+1α2e-ατb+3e-2ατhL3=α3e-2ατe-ατa23b+e-3ατa23h+3α3e-ατe-2ατa21+3e-ατa23b+7e-3ατh+1α3e-ατb+9e-2ατhwith h=aψ/2. After introducing c-=Sy, y=y1,y2,y3T, we get(82)dy1dτ=φ1τy3+d,dy2dτ=y1+φ2τy3,dy3dτ=y2+φ3τy3.The changey3=ψ(τ) reduces (82) to a linear equation of the third order:(83)ψ3-φ3τψ2-2dφ3dτ+φ2τψ1-d2φ3dτ2+dφ3dτ+φ1τψ=d.The variables y1, y2, and y3 are connected to ψ(2), ψ(1), and ψ with(84)y=TΨ,Ψ=ψ2,ψ1,ψT,T=1-φ3-dφ3dτ+φ301-φ3001.Let(85)d=ψ3-6+φ3τψ2-11+2dφ3dτ+φ2τψ1-6+d2φ3dτ2+dφ3dτ+φ1τψ.After substitution of dτ into (83), we get the equation with -1, -2, and -3 as roots of the characteristic polynomial. The return to the original variables gives (86)d=ΓτT-1τS-1τc-,Γ=-6+φ3,-11+2dφ3dτ+φ2,-6+d2φ3dτ2+dφ2dτ+φ1.It is obvious that (86) provides an exponential decrease of solutions of (80). At the final stage, we solve the IVP for the system obtained from (76) after changing the phase coordinates according to (79) with control (86). Then, we return to the original variables. The initial values for the Cauchy problem are(87)υ10=-x-1-12gx-1-124∂g∂c1x-1gx-1-g--,u20=gx-1+16∂g∂c1x-1gx-1+g-+g---,c30=0.


## 7. Numerical Simulation

For numerical simulation, we choose another problem—control of a single-link manipulator. According to [40], the system of equations describing the motion of a manipulator under perturbations has the form(88)x˙1=-x2,x˙2=-a1x2-a2sin⁡x1+u+μt,where x1 is the manipulator deviation angle in the vertical axis, x2 is the deflection angle change rate, a1=α-L-2m1-1, m1=m0+M/3, a2=gL-1m0+M/2m1-1, g is the gravity factor, α- is the viscous friction coefficient, m0 is the load mass, L is the manipulator length, and M is its mass. The vector of phase coordinates x=x1,x2T. The control u is scalar. We consider the boundary conditions(89)x0=x-,x1=0.System (13) and conditions (14) for problem (88) and (89) have the form(90)dc1dτ=αe-ατc2,dc2dτ=-αe-ατa2sin⁡c1-αe-ατa1c2+αe-ατd+μαe-ατ1-e-ατ,(91)c10=x-1,c20=0,c1τ⟶0,c2τ⟶0asτ⟶∞.To solve (90) and (91), we perform one transformation of function c2τ:(92)c2τ=c21τ-μe-ατ,which leads to new functions c2(1)(τ). The system takes the form(93)dc1dτ=αe-ατc21+αμe-2ατ,dc21dτ=-αe-ατa2sin⁡c1-αe-ατa1c21-αe-2ατa1μ+αe-ατd-αμe-2ατ.The linear part of (93) is(94)dc-dτ=αe-ατPc-+αe-ατQd,c-=c1,c21T,P=01-a2-a1,Q=01.After solving the stabilization problem (91), we obtain the control law(95)dτ=Mτc-,c-=c1,c21T,Mτ=-6αe2ατ+4α2e2ατ+2e2ατ-a2α2α2-3-3α+αe-ατa1eαταT.Next, we solve the Cauchy problem for (93) under (95) with initial data c10=x-1, c2(1)(0)=μ. Finally, we return to original variables x1, x2, u, and t. For numerical simulations, we choose x-1=1 rad, α-=0.1, α=0.25, L=10 meters, M=20 kg, m0=1 kg, and μ=0.1. Figures 1 and 2 show the control ut and the corresponding transient plots of the phase coordinates x1t, x2t, and x3t with respect to the original independent variable t.Figure 1
Time change of the phase variables.Figure 2
Time change of the controlu(t).


## 8. Numerical Simulation Discussion

Numerical simulation consists of the following stages:(1)
The auxiliary system (93) construction, performed with analytical methods and realized with computer algebra(2)
System (93) stabilization and control law (95) determination, performed with analytical methods and realized with computer algebra(3)
Change of variables fromc- to c, performed with computer algebra(4)
Solution of the IVP (90) and (91) with control (95) in c variables, performed with some numerical integration method(5)
Return to the original variables according to (8), (12), and (15) along with the IVP solution.It follows from the simulation results that(a)
control at the initial stage demands the greatest energy input which is determined by the initial stagex-1,x-2;(b)
multiple initial data simulations show that, for the given transfer time of 1 second, the reachable area of initial states is given byx-1≤π/2 and x-2≤1;(c)
the numerical simulation can be performed with an average personal computer.


## 9. Conclusions

The analysis of the theorem proof shows that the most difficult and time-consuming part of the algorithm implementation can be proceeded with analytical methods of computer algebra packages. The results of the numerical simulation of interorbital flight prove that the method can be used for construction and simulation of various technical objects control systems.



---

*Source: 8529760-2017-09-17.xml*