Computer algebra and problems of motion stability

Banshchikov Andrej V., Bourlakova Larissa A.

Institute of System Dynamics and Control Theory SB RAS

134, Lermontov St., Irkutsk, 664033, Russia

E-mail : irteg@icc.ru

Full paper in compressed Postscript *.ps.gz

The paper describes the algorithms for investigation of stability of mechanical systems. These algorithms are implemented with the aid of the computer algebra system (CAS) in the software package "STABILITY". Package investigates stability and stabilization of mechanical systems in both pure symbolic and numerical-symbolic forms [1].

Let's consider the stability task for trivial solution of matrix differential equation

$$M \ddot x\ + (2G+D)\dot x\ + (K+P)x\ =Q(x,\dot x )\ , \eqno (1)$$

where $M=M^{T}>0\ ,\ D=D^{T}\ ,\ G=-G^{T}\ ,\ K=K^{T}\ , \ P=-P^{T}\ -$are $(n\times n)$ matrices of dissipative, gyroscopic, potential and nonconservative (circulatory) forces; $\ Q(x,\dot x)\ -$ is single-column matrix of nonlinear forces, $\ Q(0,0)=0\$.

This presentation is possible for differential equations of perturbed movement in the neighborhood of stationary motion or invariant set of stationary motions. Algorithm to reduce equations to form (1) and algorithm to classify nonlinear forces $\ Q(x,\dot x )\$ suitable for implementation with CAS are described in [2]. There are three specific groups of the tasks.

Investigation of stability in the first approximation. The question of asymptotical stability or instability can be investigated by studying of the first approximation equations ( when in eq. (1) $\ Q(x,\dot x )=0$) , if we'd use classic Lyapunov's theorems.

Step 1. Construct characteristic equation of the system (1):

$$det\left(M \lambda^{2}\ + (2G+D)\lambda\ + (K+P)\ \right)=0\ . \eqno (2)$$

Coefficients of characteristic polynomial are calculated in symbolic or numerical-symbolic form. In the latter case they contain certain parameters. Values of these parameters must be selected from the condition of asymptotical stability.

Step 2. Analyze the structure of the polynomial (2).

Step 2.1. If polynomial (2) contains all powers of $\lambda$ , then we construct the Hurwitz matrix and write out all Routh-Hurwitz conditions and their analogies. It gives us necessary and sufficient conditions that all roots of (2) have negative real parts. If these conditions are met, then the motion is asymptotically stable, independently from nonlinearities $\ Q(x,\dot x )\$. These condition comprise the system of algebraic inequalities of system parameters (matrix elements of equation (1), amplification coefficients in linear controlling forces, etc). Inequalities are written out in either symbolic or numerical-symbolic form. A description of some algorithms for solving the systems of algebraic inequalities is given in [2]. The software implementation of these algorithms allows us to solve following tasks:

a). Finding the domains, which satisfy the inequalities of one-parameter polynomials with numerical coefficients

$$f_ {i} (p) = a_{n}^{(i)}\ p^n + a_{n-1}^{(i)}\ p^{n-1} + \ldo... ...\ \ \ \ (\leq 0)\ , \ \ \ (\ i = 1,\ldots, n_1\ )\ . \eqno (3)$$

The solution (the domain obtained) is represented as a union of disjoint intervals of the numerical axis.

b). Search and graphical constructing of two-dimensional domains corresponding to a system of algebraic inequalities

$$f_{i}\ (p_1\ ,\ p_2) \ \geq 0 \ \ \ \ \ (\leq 0) \ , \ \ \ (\ i = 1,\ldots,n_2\ ) \ , \eqno (4)$$

where $\ p_1\$ is the parameter reckoned on the abscissa axis, and $\ p_2\$ is the parameter reckoned on the axis of ordinates. The boundaries of variation of the parameters $\ p_1\$ and $\ p_2\$for constructing 2D - domains are determined as the solutions of problem (a).

Step 2.2. If in polynomial (2) some powers of $\lambda$ are absent, we must conclude that the system is structurally unstable.

Step 2.2.1. Let in polynomial (2) all coefficients be $\ a_{2k-1}=0\ \ (k=1,\ldots , n) .$This is always true if in linear part of the system are present (i) positionary forces (potential or non-conservative) or (ii) together with potential forces are present gyroscopic ones. Stability of such a system is possible only in a critical case when all roots of characteristic polynomial are imaginary or(and) zero. Let's substitute $\lambda ^{2}=\nu \$ in equation (2). For resulting equation of power n let's write out the conditions when all roots of the polynomial are negative (different) real numbers [3]. This way we get the stability conditions of linear system and necessary conditions of stability of nonlinear system (1).

Direct calculation of characteristic equation (2) requires us to calculate large determinants. Thus it's preferable to use some methods of decomposition. Let's suppose that lagrangian $\ L\$ of the system allows presentation in the form

$$L(x_{1},\ldots , x_{n},\dot {x}_{1},\ldots , \dot {x}_{n}) = \tilde {L} \ +\ \tilde{\tilde{L}}\ ,$$

$${\rm where} \ \ \ \tilde{L}= \tilde{ L}(x_{1},\ldots , x_{m},... ...\ldots , \dot {x}_{m},\dot {x}_{m+1},\ldots , \dot {x}_{k})\ ,$$

$$\tilde{ \tilde{L}}= \tilde{\tilde{ L}}(x_{m+1},\ldots , x_{k}... ...\dot {x}_{k},\dot {x}_{k+1},\ldots , \dot {x}_{n})\ . \eqno(5)$$

For the lagrangian (5) let's construct Lagrange equations of second kind, that can be transformed to form (1). When constructing characteristic equation for the system, where lagrangian is quadratic with respect to coordinates and velocities, we must calculate determinants of dimension $\ (k-m),\ (m), \ (n-k)\ .$ Indeed, considering (5), characteristic matrix $\ B\$ of the system has blocked structure:

$$B=\left (\begin{array}{ccccc} k-m & m & n-k & \\ \overbrac... ... \\ B_{31} & 0 & B_{33} & \} & n-k \end{array} \right )\ .$$

If matrices $\ B_{22}\$ and $\ B_{33}\$ are nondegenerate, then

$$det B=\left \vert\ B_{11}\ -\ B_{13}B^{-1}_{33}B_{31}\ - \ B_... ...left \vert B_{22}\right \vert \left \vert B_{33}\right \vert .$$

This formula can be generalized for the case when the matrix form of characteristic equation has the form

$$\left \vert\begin{array}{ccccc} B_{11} & B_{12} & B_{13} & \l... ...\ldots & B_{NN} \end{array} \right \vert \ =\ 0 .\ \ \eqno(6)$$

Here $\ B_{ii}\ (i=1,\ldots , \ N)\$ are $\ (n^{i}\times n^{i})\$matrices $\ ( \ \sum n_{i}=n\ ) .$For example, this is the structure of characteristic polynomial for differential equations, describing a motion of system of connected bodies in potential field. If $\ det B_{ii}\neq 0\ \ \ (i=2,\ldots , N)\ ,$ then (6) can be written out in the form

$$\left \vert\ B_{11}\ -\ B_{1N}B^{-1}_{NN}B_{N1}\ - \ldots \ -... ... B_{22}\right \vert\ldots \left \vert B_{NN}\right \vert =0\ .$$

Consequentially, calculation of characteristic determinant is reduced to calculation of several determinants of significantly lesser dimension. It allows us to study the systems with much more degrees of freedom.

Investigation of stability by second Lyapunov's method. Sufficient conditions of stability, asymptotic stability or instability can be obtained by using theorems of second Lyapunov's method.

In critical cases , only the Lyapunov's second method is used. But then there arises the problem of the proximity of sufficient stability conditions to necessary conditions, i.e. of a "good" choice of the Lyapunov function . As investigations show, in critical cases sufficiently soft stability conditions are obtained when the first integral of the differential equations of motion is taken as the Lyapunov function.

Definition. The relation $\ df(\dot {x}, x)/dt=F_ {1}(\dot {x}, x)\$is said to be the differential consequence of equations (1) if it holds for all solutions of the system.

In the theorems of the Lyapunov's second method, we shall use differential consequences of equations (1). One of the algorithms for constructing differential consequences is described in [1] and involves the following procedure : we multiply (1) on the left by matrices of the form $\ x^ {T} L_ {1}\ , \ \dot x^ {T} L_ {2}\ , \ \ddot x^ {T} L_ {3}\ ,$ and transform the final expressions by selecting the symmetric parts. The matrices $\ L_ {i} \$ are equal to those of the initial equation or their combinations.

Example 1. Let $\ L_ {3} =GM^ {-1}\ ,\ L_ {2} =KM^ {-1}\$, and $\ L_{1} =GM^ {-1}\$. After selecting the full derivative, we obtain three consequences. We combine them into one bunch :

$${ d\over dt} (\dot {x} ^ {T} (4G^ {T} M^ {-1} G + K) \dot {x}... ...-1} G) x + 2 {x} ^ {T} (KM^ {-1} G + PM^ {-1} G-G) \dot {x}) =$$

$$= -2 \dot {x} ^ {T} (GM^ {-1} G^ {T} M^ {-1} D + KM^ {-1} D +... ...P) \dot {x} -2x^ {T} (G^ {T} M^ {-1} K + G^ {T} M^ {-1} P) x +$$

$$+ 2x^ {T} (2KM^ {-1} G^ {T} M^ {-1} D - 2PM^ {-1} GM^ {-1} D-G^ {T} M^ {-1} D-P^ {T} M^ {-1} K) \dot x\ . \eqno (7)$$

Now, by $\ W\$ we denote the form in (7) on the right. If in (7) the form under the sign of derivative has the fixed sign and the form $\ W\$ has the fixed sign with an opposite sign, in accordance with Lyapunov's theorem the trivial solution of the system is asymptotically stable irrespective of nonlinear forces. If the form $\ W\$ in (7) is of constant signs (the sign opposite to that of the form on the left), and the set $\ {\cal L} =\ \{\ x,\ \dot x\ :\ W=0\ \}$does not contain the full trajectories except for $\ x=\dot x=0\$, then the trivial solution of the linear system (1) is asymptotically stable. $\lhd$

Let's define the function $\ V\$ :

$$V=y^{T}{\cal N}y + aF(y)= \dot{x}^{T}N \dot{x} +x^{T}Lx+x^{T}B^{T}\dot{x}+ \dot{x}^{T}Bx+aF(x, \dot{x})\ , \eqno(8)$$

here

$$L=L^{T}\ ,\ \ \ N=N^{T}\ ,\ \ \ \ {\cal N}= {\cal N^{T}}= \le... ...= \left ( \begin{array}{c} x \\ \dot x \end{array}\right )\ ,$$

a is constant that can be equal to 0 or 1; $\ F(y)$ - function of power bigger than two.

Let's calculate derivative of function (8) for equation (1):

$$\dot V =y^{T}( A^{T}{\cal N }\ + \ {\cal N }A)y \ + \tilde{Q}... ...\ + a \frac{\partial F}{\partial y}(Ay\ +\ M^{-1}\tilde{Q})\ .$$

here

$$A= \left ( \begin{array}{cc} 0 & E \\ -M^{-1} (K+P) & -M^{-1... ...de{Q}= \left ( \begin{array}{c} 0 \\ Q \end{array}\right )\ .$$

We require that this derivative was presented in the form:

$$\dot V\ =\ \alpha V\ +\ y^{T} Wy\ + \Phi (y)\ .$$

Here $\ \Phi (y)$ is nonlinear function of power bigger than two. Then we obtain the equation system:

$$A^{T}{\cal N }\ + \ {\cal N }A \ -\ \alpha {\cal N}\ = W\ ; \eqno(9)$$

$$Q^{T}M^{-1} ( N\dot x + Bx ) + ( \dot x^{T}N + x^{T}B^{T})M^{-1}Q \ +$$

$$+\ a\ (\ \frac{\partial F}{\partial x}\dot x \ - \ \frac {\pa... ... \ (K + P)x \ - \ Q )\ )\ = \ a\alpha F\ +\ \Phi (x,\dot x)\ .$$

From last equation we find $\Phi (x,\dot x)\ .$For equation (9) we will consider several possibilities.

Variant 1. $\ W\$ is given definitively positive (negative) matrix. From classical theorems it is known that if at least one of characteristic equation (2) roots has positive real part, for equation (9) there always exist such $\ \alpha >0 \$ and $\ {\cal N}$ that $\ y^{T}{\cal N}y \$ isn't sign-definite with sign, opposite to sign of $\ W\ .$ Vice versa, if for equation (9) with $\ W>0\$ are found $\ \alpha >0 \$ and $\ {\cal N}$ and $\ y^{T}{\cal N}y \$ is not definitively negative, then trivial solution of equation (1) is not stable and some roots $\lambda$ of eq. (2) have positive real part.

Variant 2. $\ \alpha =0\ ,\ W\$ is given definitively negative (positive) matrix. It is known that if all roots of characteristic equation (2) of system (1) have negative real parts, then in (9) there is only one solution for $\ {\cal N}$. In these conditions obtained matrix is definitively positive (negative). Vice versa, if with given $\ W<0 \ (>0)\$ it's found solution $\ {\cal N}$ of equation (9) and $\ y^{T}{\cal N}y>0\ (<0)\$, then zero solution of equation (1) is asymptotically stable and all $\lambda$ of equation (2) have negative real parts.

If roots $\lambda$ of equation (2) are such that $\ \sum m_{i} \lambda _{i} \neq 0\ \ (\ i=1,\ldots , n,\ \sum m_{i}=2\ )$ and some of them have positive real part, then equation (9) have single solution $\ {\cal N}$ and $\ y^{T}{\cal N}y \$ is not definitively positive (negative). Zero solution of system (1) is unstable.

Variant 3. $\ W=0\ .$ Then equation (9) is equivalent to equation [4] :

$${\cal C}{\it x}= 0\ , \eqno(10)$$

where

$${\it x}= \left ( \begin{array}{c} {\cal N}_{1*}^{T} \\ . \\... ...)\ ,\ \ \ \ {\cal C}= (A-\alpha E)\otimes E\ +\ E\otimes A \ .$$

This equation has nontrivial solution $\ \Longleftrightarrow \alpha =\lambda _{r}+ \lambda _{s}\ ,$ where $\ \lambda_{r},\ \lambda _{s}\$ are eigenvalues of matrix $\ A\ .$

Let's write out equation (10) in the form :

$$B\ +\ B^{T}\ -\ (D-2G)M^{-1}N\ -\ NM^{-1}(D+2G) \ -\ \alpha N\ =\ 0\ ;$$

$$(K-P)M^{-1}N - NM^{-1}(K+P) + B^{T}M^{-1}(D+2G) - (D-2G)M^{-1}B + \alpha (B^{T}-B) = 0\ ;$$

$$\alpha ((K-P)M^{-1}N\ +\ NM^{-1}(K+P)\ +\ B^{T}M^{-1}(D+2G)\ + \ (D-2G)M^{-1}B)\ +$$

$$+\ \alpha ^{2}(B^{T}+B)\ +\ 2((K-P)M^{-1}B\ +\ B^{T}M^{-1}(K+P))=\ 0\ ;$$

$$2L = (K-P)M^{-1}N + NM^{-1}(K+P)+ B^{T}M^{-1}(D+2G) + (D-2G)M^{-1}B + \alpha (B+B^{T})\ . \eqno(11)$$

If eigenvalues of matrix $\ A$ are known and none of them have real part equal to zero, then stability task can be solved without using Lyapunov functions. Let's put limitations on structure of $\ V\ .$

Let $\ B=0\ .$ Then matrix $\ N$ is a solution of over-defined system :

$$(D-2G)M^{-1}N\ +\ NM^{-1}(D+2G) \ +\ \alpha N\ =\ 0\ ;$$

$$(K-P)M^{-1}N\ -\ NM^{-1}(K+P)\ =\ 0\ ;$$

$$(K-P)M^{-1}N\ +\ NM^{-1}(K+P)\ =\ 0\ . \eqno(12)$$

If system (12) has a solution $\ N^{*}\$, then $(N^{*}+N^{*T})\$, $\ (N^{*}-N^{*T})\$and $\ ({\cal Q}^{T}N^{*}{\cal P}\ + \ {\cal P}^{T}N^{*}{\cal Q})\ ,$( where $\ {\cal P},\ {\cal Q}\$ are commutative matrices of $\ M^{-1}(D+2G)\ ,\ M^{-1}(K+P)\$) are solutions too.

From first equation it follows that $\ \alpha=-(\mu_{r}+\mu_{s})\ , \ \mu_{r},\ \mu_{s}\$ are eigenvalues of matrix $\ M^{-1}(D+2G)\ .$ Second equation always have set of nontrivial solutions. From second and third equations if follows the condition $\ det (K+P)= 0\$ when system (1) is critical. Then third equation (12) has nontrivial solution. The solution of system (12) can be used as Lyapunov function $\ V=V(\dot x)\ ,$ because if $\ N$ satisfies system (12) then $\ L=0\ .$

Example 2. Consider a system with two degrees of freedom :

$$\ddot x_{1}+ c_{11}\dot x_{1}+ c_{12}\dot x_{2}+ q_{11}x_{1}+q_{12}x_{2}=Q_{1}\ ;$$

$$\ddot x_{2}+ c_{12}\dot x_{1}+ c_{22}\dot x_{2}+ q_{12}x_{1}+q_{22}x_{2}=Q_{2}\ ;$$

Eigenvalues of matrix $\ (D+2G)=C\$ are solutions of the equation

$\ \mu^{2}-\mu (c_{11}+c_{22})+c_{11}c_{22}-c_{12}c_{21}=0\ .$System (12) has the form :

$$n_{11}(\alpha + 2c_{11})+n_{12}c_{21}=0\ ; \ \ \ n_{22}(\alpha + 2c_{22})+n_{12}c_{12}=0\ ;$$

$$n_{11}c_{12}+n_{12}(\alpha + c_{11}+ c_{22})+n_{22}c_{21}=0\ ;$$

$$n_{11}q_{11}+n_{12}q_{21}=0\ ;\ \ \ \ n_{12}q_{12}+n_{22}q_{22}=0\ ;$$

$$n_{12}(q_{11}+q_{22})+n_{11}q_{12}+n_{22}q_{21}=0\ ;\ \ \ \ \ n_{12}(q_{11}-q_{22})-n_{11}q_{12}+n_{22}q_{21}=0\ ; \eqno(13)$$

Value for $\alpha$ must be chosen from $\ -(\mu _{r}+\mu_{s}).$

Solution 1. Let $\ \alpha = -(c_{11}+c_{22})\ .$ When $\ c_{12}=-c_{21}=2g\ , \ c_{11}=c_{22}=c\ ,\ q_{11}=q_{22}=q_{12}=q_{21}=0\$ system (13) allows the solution : $\ n_{11}=n_{22}=n\ ,\ n_{12}=0\ .$ Consequently (8) has the form :

$$V(\dot x)=n(\dot x_{1}^{2}+ \dot x_{2}^{2})>>0 \ \ (n>0)\ .$$

$$\dot V= V_{1}=-2cV+2n(Q_{1}\dot x_{1}+ Q_{2}\dot x_{2})\ .$$

If system (13) contains only dissipative and gyroscopic forces, then for $\ c > 0\$ trivial solution is asymptotically stable with respect to velocities despite any nonlinear forces $\ Q=Q(\dot x)$ because $\ V_{1}<<0\ .$If $c<0\$ trivial solution is unstable because $\ V_{1}>0\ .$

Solution 2. Let $\ \alpha = -2c_{11}\ .$ Then from system (13) we obtain : $\ c_{12}=d_{12}+2g=0\ ,\ q_{11}=k_{11}=0\ ,\ q_{12}=k_{12}+p=0\$, $\ q_{22}\neq 0\$, $\ q_{21}=k_{12}-p\neq 0\ ,\ c_{21}=d_{12}-2g\neq 0\ , \ c_{11}\neq 0\$( in the system are present dissipative, gyroscopic, potential and nonconservative positional forces), $\ n_{12}= 0\ ,\ n_{22}=0\ ,\ n_{11}\neq 0\ .$Lyapunov function has the form $\ V=\dot x_{1}^{2}\$. For eqs. (13) $\ dV/dt= V_{1}=-2c_{11}\dot x^{2}+ \dot x_{1}Q_{1}\ .$ When $\ c_{11}>0 \$and $\ Q_{1}= Q_{1}(\dot x_{1})\$ trivial solution is asymptotically stable with respect to $\dot x_{1}\ .$ When $\ c_{11}<0 \$ trivial solution is unstable. $\lhd$

Variant 4. $\ \alpha =0\ ,\ W=0\ .$ Then we obtain two groups of conditions when equation (1) has first integral of form (8) :

a)

$$B\ +\ B^{T}\ -\ (D-2G)M^{-1}N\ -\ NM^{-1}(D+2G) \ =\ 0\ ;\eqno(14)$$

$$(K-P)M^{-1}N\ -\ NM^{-1}(K+P)\ +\ B^{T}M^{-1}(D+2G)\ -\ (D-2G)M^{-1}B\ =\ 0\ ; \eqno(15)$$

$$(K-P)M^{-1}B\ +\ B^{T}M^{-1}(K+P)=\ 0\ ; \eqno(16)$$

$$2L\ =\ (K-P)M^{-1}N\ +\ NM^{-1}(K+P)+\ B^{T}M^{-1}(D+2G)\ +\ (D-2G)M^{-1}B\ .$$

b)

$$\frac{\partial F}{\partial x} +2NM^{-1}Q-(D-2G)M^{-1}\frac{\partial F}{\partial \dot{x}}=0\ ;$$

$$2B^{T}M^{-1}Q-(K-P)M^{-1}\frac{\partial F}{\partial \dot{x}}=... ... \ \ \ \ \ Q^{T}M^{-1}\frac{\partial F}{\partial \dot{x}}=0\ .$$

Equations of group (a) give necessary and sufficient conditions of first integral existence fir linear system (1) and existence of quadratic part of the integral of form (8) for nonlinear system.

Equations of group (b) are used to build $F(x,\dot{x})$ if it exists or to find such $\ Q\$ that eq. (8) is an integral. Function $F(x,\dot{x})$ must fulfill the conditions:

$$\frac{\partial^{2} F}{\partial \dot{x}_{i}\partial \dot{x}_{j... ...al x}= \frac{\partial ^{2} F}{\partial x\partial \dot{x}}\ \ .$$

$$(\ i,\ j\ = 1,\ldots ,n\ )$$

System (a) (14)-(16) always has trivial solution $\ N=0,\ B=0\ (\ B^{T}=0\ )$. It is a system of uniform linear equations and can have (i) only trivial solution or (ii) set of solutions. If there exists nontrivial solution $\ N^{*}\ ,\ B^{*}\$ then $\ N_{1}^{*}=N^{*}\pm N^{*T}\ ,\ B_{1}^{*}=B^{*}\pm B^{*T}\$, $\ N^{**}={\cal P}N^{*}{\cal Q}^{T}+{\cal Q}N^{*}{\cal P}^{T}\ ,\ B^{**}={\cal P}^{T}B^{*}{\cal Q}+{\cal Q}^{T}B^{*}{\cal P}\$ ( where ${\cal P},\ {\cal Q}\ -$ are matrices commutative with $\ M^{-1}(D+2G)\ ,\ M^{-1}(K+P)\$ ) are solutions too. Existence conditions for nontrivial solution give us limitations on matrices of original system.

Example 3. Let $\ D=0\ ,\ P=0\ ,\ Q=0\$ ( system with potential and gyroscopic forces). Solution $\ N=K-4GM^{-1}G\ ,\ B=2G^{T}M^{-1}K\$ determines first integral

$$\ V_{1}=\dot{x}^{T}K\dot{x} +(Kx +2G\dot{x})^{T}M^{-1}(Kx+2G\dot{x})= const\ . \ \ \ \ \lhd \eqno(17)$$

Let's outline some properties of matrix equations of special form that can be obtained by transformation of equations (14)-(16) :

(i) equation $\ AX-XA^{T}=0 \ \ (\ A,\ X\$ are ${n\times n}\$ matrices ) always has symmetric nontrivial solution; skewed-symmetric solution exists if matrix $\ A$ has multiple eigenvalues with corresponding prime elementary divisors and/or different Jordan blocks;

(ii) equation $\ AX+XA^{T}=0\$ has nontrivial solution if some of roots of characteristic polynomial $f(\mu )\$ of matrix $\ A$ are equal to zero and/or $f(\mu )\$ is even function.

Let's discuss obvious variants for system (a) (14)-(16).

Variants 4.1. $\ B=-B^{T} \ .$ The solution of $\ N$ would be nontrivial if characteristic polynomial of matrix $\ M^{-1}(D+2G)\$ must be even or $\ det(D+2G)=0\$; characteristic polynomial of matrix $\ M^{-1}(K+P)\$ must have multiple roots ( if $\ B\neq 0 \$). Consequence of eqs. (14)-(16) is the equation:

$$(K-P-D+2G)M^{-1}(N+B) - (N+B)M^{-1}(K+P+D+2G) = 0\ , \ \ (\ A_{1}X=XA_{2}\ ) \eqno(18)$$

Equation (18) has nontrivial solution if some of eigenvalues $\chi\$ of matrix

${\cal C}_{1}=(K-P-D+2G)M^{-1}\otimes E\ -\ E\otimes (K-P+D-2G)M^{-1}\$are equal to zero; $\ \chi = \mu _{s}-\nu _{r}\ ,$ where $\ \mu_{s}$ are eigenvalues of matrix $\ M^{-1}(K+P-D-2G)\$, $\ \nu _{r}\$ are eigenvalues of matrix $\ M^{-1}(K+P+D+2G)\ .$ The number of solutions is $\ n_{1}=\sum \delta _{sr}\$, where $\ \delta _{sr}\$ is a power of maximal common denominator $\ (\mu -\mu_{s})^{p_{s}}\$ and $\ (\mu -\nu_{r})^{q_{r}}\$of characteristic polynomials of these matrices. Let's build equations $\ det{\cal C}_{1}\ = \ 0\$ and find existence conditions for $\chi\$equal to zero. If the equation (18) has a solution $\ X^{*}\$ then also exist a solution $\ X^{**}={\cal P}_{1}X^{*}{\cal P}_{2}\$, where matrix $\ {\cal P}_{1}\$ is commutative with $\ A_{1}\$, $\ {\cal P}_{2}\$ is commutative with $\ A_{2}\ .$

Let's select symmetrical and skewed-symmetrical parts of solution $\ X=N+B\$ and substitute them to equations (14) ¨ (16). Then we obtain conditions on system matrices.

Example 4.

4.1. Let $\ K=0\ ,\ P=0\ ,\ D=0\ ,\ Q=0\$ (system with gyroscopic forces). System (14)-(16) has solution $\ B=G\ ,\ N=M\$. Equation (1) has first integral $\ \dot x^{T}M\dot x\ -\ 2x^{T}GM^{-1}Gx\ +\ 2\dot xGx=const.$

4.2. Let $\ K+P=M{\cal K}\ ,\ (\ {\cal K}={\cal K}^{T}\ )\ ,\ \ D+2G=M\Gamma\ , \ (\ \Gamma =-\Gamma ^{T}\ ).$ Then matrices $\ A_{1},\ A_{2}\ (18)\$ are identical. The task is reduced to calculation of matrices commutative with $\ (\ {\cal K}+\Gamma\ ).$ Also, $\ {\cal K}B=B{\cal K}\ , \ \Gamma N=N\Gamma\ .$Note that if $\ {\cal K}$ does not have multiple eigenvalues then $\ B=0\ .$$\lhd$

Variant 4.1.1. Let $\ N=0\ .$ Then from eq. (14) follows $\ B=-B^{T} \ .$This solution is possible if above-mentioned conditions are true: characteristic polynomial of matrix $\ M^{-1}(D+2G)\$ must be even or $\ det(D+2G)=0\$; characteristic polynomial of matrix $\ M^{-1}(K+P)\$ must have multiple roots.

Example 5.

5.1. Let $\ D=0\ ,\ P=0\ ,\ K=\alpha M\ ,\ Q=0\$ ( system with gyroscopic and potential forces ). System (14)-(16) has solution $\ B=G ,\ N=0\$That gives first integral of the equation (1) $\ V_{3}= 2\dot x^{T}Gx\ +\ x^{T}GM^{-1}Gx= const\ .$

5.2. Let $\ K=P=0\ , \ D\neq 0\ ,\ G\neq 0\ ,\ Q=0\ .$Solution $\ B=G ,\ N=0\$ with $\ GM^{-1}D+DM^{-1}G=0\$ determines first integral of equation (1): $\ V_{4}= 2x^{T}G\dot x\ +\ x^{T}GM^{-1}(D+2G)x=const.$ $\lhd$

Variant 4.2. Let $\ B=0\$. Then it remains a system of two equations :

$$(D-2G)M^{-1}N + NM^{-1}(D+2G) = 0\ ; \ \ \ (K-P)M^{-1}N - NM^{-1}(K+P) =0\ . \eqno(19)$$

First equation has nontrivial solution if $\ det(D+2G)=0\$ or characteristic polynomial of that matrix is even function. Second equation always has symmetric solution. The solution of this system satisfies the equation $\ (K-P-D+2G)M^{-1}N - N(K+P+D+2G) = 0\$ and existence conditions for nontrivial solution, similarly to eq. (18).

Example 6.

6.1. Let $\ D=0\ ,\ P=0\ ,\ Q=0\$ ( system with potential and gyroscopic forces ). Then in eq. (18) $\ A_{1}=A^{T}_{2}\$, characteristic polynomials of matrices are identical. There exist a solution $\ N=M\ ,\ B=0\$ that satisfies eqs. (14) and (16). The system allows a first integral $\ V_{2}=\dot x^{T}M\dot x \ +\ x^{T}Kx\ = const.$From theorems of Lyapunov's second method it follows that if $\ K>0\$then trivial solution of linear system (1) is stable with respect to variables $\ \dot x, x\$. If matrix K is not definitely positive, we generate a bunch of forms $\ V_{2}\$ and $\ V_{1}\ (17)\$. Conditions of sign-definitiveness of this band give sufficient conditions of stability for the critical system in question.

6.2. Let $\ K=P=0\ , \ Q=0\$( system with dissipative and gyroscopic forces). System (19) has solution $\ N=GM^{-1}D\ +\ 2GM^{-1}G\$ where

$$\ GM^{-1}D+DM^{-1}G=0\ .\eqno(20)$$

First integral of equation (1) has the form: $\ V_{5}= \dot{x}^{T}(GM^{-1}D\ +\ 2GM^{-1}G)\dot x=const\ .$If $\ V_{5}<<0\$ then trivial solution of linear equation (1) is stable with respect to $\dot x\ .$Note that condition (20) is met only when $\ D\$ is matrix of variable sign. In particular, if $\ M=E\$ then $\ trD=0\ .$

6.3. Let $\ P=0\ , \ Q=0\ ,\ K=\beta K_{1}\ ,$ where $\ K_{1} = GM^{-1}D + 2GM^{-1}G\$(system with dissipative, gyroscopic and potential forces) and condition (20) is met. Linear equation (1) has first integral

$$V_{6}= \dot{x}^{T}K_{1}\dot x\ +\ \beta x^{T}K_{1}M^{-1}K_{1}x\ =const\ .$$

If $\ (GM^{-1}D\ +\ 2GM^{-1}G)<<0\ ,\ \beta <0\$ then $\ V_{6}<<0\ ,$ and trivial solution of eq. (1) is stable with respect to $\ x,\ \dot x.$$\lhd$

Other variants of solutions of systems (a) and (b) are discussed in [2]. Integral of the linear system can be used as Lyapunov function for nonlinear system too.

The calculated first integrals (or their bunches) can be used to construct Lyapunov functions in proving the stability of the trivial solution for the system (1) if Sylvester's inequalities of sing-definiteness of the quadratic part of the integral (8) with respect to all variables of the problem are satisfied.

Systems with decomposition. Let system (1) be separable to primary and secondary subsystems:

$$\left( \begin{array}{cc} M_{11} & M_{12} \\ M_{12}^{T} & M_{... ...egin{array}{c} \dot x_{1} \\ \dot x_{2} \end{array}\right) \ +$$

$$+\left( \begin{array}{cc} P_{11} & P_{12} \\ -P_{12}^{T}& P_{... ...\begin{array}{c}Q_{1} \\ Q_{2} \end{array} \right) \ \eqno(21)$$

For such representation of a system we can write out differential consequences and use them in more "subtle" theorems than Lyapunov theorems on stability and asymptotic stability.

For example, if $\ P=0 \$ then from a differential consequence

$${d\over dt}(\dot{x}^{T}M\dot{x} + x^{T}Kx ) = 2\dot{x}^{T}Q - 2 \dot{x}^{T}D\dot{x} - 2\dot{x}^{T}Px \$$

by Krasovsky theorem it follows a modification of theorem in [5] :

Theorem . If $\ D_{11}>0\ , \ D_{22}\geq 0\ , \ D_{12}=0\ , \ K>0\$and set ${\cal L} = \{(x_{1},x_{2},\dot x_{1},\dot x_{2})\ :\ \dot x_{1}^{T}D_{11}\dot x_{1}+ \dot x_{2}^{T}D_{22}\dot x_{2}=0\}\$does not contain full non-zero trajectory, then zero solution of a system (21) (P=0) is asymptotically stable with respect to $\ \dot x, x\$ with any non-linear gyroscopic and potential forces. $\Box$

If $\ Q=0\$ then conditions of non-existence of full trajectory on a set $\ {\cal L}\$ lead to the statement that a system

$$M_{12}\ddot x_{2} -2G_{12}\dot x_{2}+K_{12}x_{2}=0$$

$$M_{22}\ddot x_{2} + (D_{22}+2G_{22})\dot x_{2}+K_{22}x_{2}=0$$

has only a zero solution $\ x_{2}=0\$ [5].

As it can be seen from presented theorems, studying of a system stability is reduced to searching of the conditions of sign-definitiveness, sign-constantness and equality to zero for matrices. Note that the task of sign-constantness for matrices cannot be solved on a computer without symbolic calculations. This is why existing numerical packages for stability investigation can work only with rough systems. This is also true for application of N.N. Krasovsky and other theorems based on Lyapunov functions with constant sign derivatives, where we must solve a question about absence of full trajectories of the system on some manifold.

Let's present one of possible algorithms for solution of this task.

1). Choosing a Lyapunov function $\ V(x,\dot x)\$ and calculation $\ \dot V\$ by equations of disturbed motion (1).

$1^{\prime}).\$ Choosing a differential consequence .

2). Checking sign-constantness of a function $\ \dot V\$ ( if $\ \dot V \geq 0\ ,$ then a group of values $(x,\ \dot x)$ when $\ \dot V = 0\$ will define a manifold we seek for).

3). If a manifold represents some surface

$F(x_1, \ldots,\ x_n,\ \dot x_1, \ldots ,\ \dot x_n) = 0\ ,\$then matrix condition:

$${\partial F \over \partial x}\dot x + {\partial F \over \part... ...dot x} (\ M^{-1}Q - M^{-1}(D+2G)\dot x - M^{-1}(K+P)x\ )\neq 0$$

(i.e. if left part isn't equal to zero at least sometimes) is a sufficient for non-existence of full trajectories in the manifold.

Note :

a). if during solution of specific tasks a sign-constant function $\ \dot V\$ can be factorized to $\ \dot V = (\ z_1\ )^{m_1}(\ z_2 \ )^{m_2}\ldots (\ z_k\ )^{m_k}\ ,\$then instead of function $\ F\$ (which determines a manifold equation) we must use all independent situations $\ \ F_i = z_i = 0\ , \hskip .4cm (\ i=1,\ldots, j\ ,\ \ j \leq k \ )$and investigate a question about non-existence of full trajectories for each of them.

b). if $\ \dot V\$ isn't a function of all variables in the task, like

$\ \dot V = \dot V(x_1,\ldots,\ x_{i_1},\ \dot x_1,\ldots,\ \dot x_{i_2})\ ,\$ £¤¥ $\ i_1+i_2 < 2n\$ then we can use $\ \ x_1=\ldots =\ x_{i_1}=\ \dot x_1=\ldots =\ \dot x_{i_2} = 0\ ,$ $\ x_{i_1+j}\neq 0 \ \ (j=1,\ldots,n-i_1)\ , \ \ \dot x_{i_2+j}\neq 0\ \ \ (j=1,\ldots,n-i_2)\$ as the manifold for investigation.

This algorithm is implemented in programming language of "Mathematica" system [6].

An example of using the "STABILITY" package. We conduct parametric analysis of asymptotic stability conditions of the equilibrium position of a satellite with controlled gravitational stabilizer on circular orbit [7] (not considering elasticity of the rod and without a gyrowheel in a satellite). Linear equations of motion are subdivided to two subsystems. In pitch angle ( $\ \theta \$) they have form

$$c\ \ddot \theta + f\ \ddot \alpha + 3\omega^2\ [\ (b-a)\ \theta +f\ \alpha\ ]\ =\ 0$$

$$f\ \ddot \theta + g\ \ddot \alpha + 3\omega^2 f\ (\ \theta +\alpha\ )\ =\ Q_{\alpha}\ \ , \eqno(22)$$

where $\ \omega\$ is an absolute value of orbital angular velocity;

$$b = J_x + m\ r\ (r+d) + \frac{1}{3}\ m\ d^{\ 2} + m_0\ (r+d)^... ... J_z + m\ r\ (r+d) + \frac{1}{3}\ m\ d^{\ 2} + m_0\ (r+d)^2\ ;$$

$$a = J_y\ ;\ \ \ g = (\frac{1}{3}\ m + m_0\ )\ d^{\ 2}\ ; \ \ ... ...ac{1}{2}\ m\ r\ d + \frac{1}{3}\ m\ d^{\ 2} + m_0\ d\ (r+d)\ ;$$

$m\ ,\ m_0\$ are, respectively, masses of the rod and the weight on its end ; $\ d > 0\$ is the length of the rod ; $\ r > 0\$ is the distance between the mass center of the system and attachment point of the rod; $\ J_x\ ,\ J_y\ ,\ J_z\$ are the principal moments of inertia of the satellite ; $\alpha$ is the rotation angle of the rod relative to the satellite body.

If $\ Q_{\alpha} = Q_{\beta} = 0\$, then potential forces influence the system. They are determined by forces of gravitational attraction and by the motion along the orbit. Such a system cannot be asymptotically stable. Indeed, the differential equations allow for the integral of energy (example 6.1) :

$$c\dot \theta ^{2}+g\dot \alpha ^{2}+2f\dot \theta \dot \alpha +3\omega ^{2}(b-a)\theta ^{2}+ 3\omega ^{2} f\alpha ^{2} =const.$$

The requirements of definite positiveness of this quadratic form give conditions for secular stability of the system.

Let's consider following control actions at the attachment point:

$$\ Q_{\alpha} = k_{\theta}^{\prime}\ \dot \theta - k_{\alpha}^{\prime}\ \dot \alpha\ .$$

The system is being influenced by potential forces (with matrix $\ K\$), gyroscopic forces (with matrix $\ G\$) and dissipative forces (with matrix $\ D\$) :

$$M\ =\ \left( \begin{array}{cc} c & f \\ f & g \end{array} \ri... ...eft( \begin{array}{cc} b-a & f \\ f & f \end{array} \right)\ ;$$

$$G\ =\ \left( \begin{array}{cc} 0 & {k_{\theta}^{\prime}}/2 \\... ...heta}^{\prime}}/2 & k_{\alpha}^{\prime} \end{array} \right)\ .$$

Let parameters of the system be : $\ \ d = 12$ m ; $\ \ r = 2.5$ m ; $m_0 = 1.531\ kg\ ;$

$m = 0.0143\ kg\ ;\ \ \omega = 0.0011$ rad/s. ; $\ \ J_y = 28\ kgm^2\ ;$ $J_x = 75\ kgm^2\ ; \ \ J_z = 76\ kgm^2\ .$With these values, domain where parameters $\ k_{\theta}^{\prime}\ ,\ k_{\alpha}^{\prime}\$ could be chosen is presented on figure 1.

Figure 1.

This domain is found by program $\ Region\$ that plots 2D solutions of numerical-symbolic algebraic inequalities (in this case Routh-Hurwitz conditions).

It is known that the problem about damping of transient process in the system is connected with evaluation of the degree of stability (i.e. with the evaluation of a disposition of the roots of characteristic equation in left-hand half-plane). Also there is an important task of choosing the parameters of the system in question, so its degree of stability will be larger or equal than a given value $\ \varepsilon\$. Let's study this question for the system (22). For the characteristic equation

$$det\ \Vert\ M\ (\ \mu - \varepsilon\ )^2 + D\ (\ \mu - \varep... ...,\ k_{\theta}^{\prime},\ k_{\alpha}^{\prime}) \ \mu^i\ =\ 0\ ,$$

( where $\ \varepsilon > 0\ -\$ is the degree of stability (offset from imaginary axis)) with the help of package "Stability" we obtained Routh-Hurwitz conditions and made parametric analysis of these conditions. The analysis showed, that the domain for choosing parameters $\ k_{\theta}^{\prime}\ ,\ k_{\alpha}^{\prime}\$exists only when $\ \ \varepsilon\ \in \ \ ]\ 0\ ,\ 0.00035\ ]\$. For example, when $\ \varepsilon = 0.0001\$

$${\tilde a}_0\ (\ k_{\alpha}^{\prime},\ k_{\theta}^{\prime}\ )... ... 1.\ k_{\alpha}^{\prime} - 0.710183\ k_{\theta}^{\prime}\ )\ ;$$

$${\tilde a}_1\ (\ k_{\alpha}^{\prime},\ k_{\theta}^{\prime}\ )... ...+ 1.\ k_{\alpha}^{\prime} + 0.69747\ k_{\theta}^{\prime}\ )\ ;$$

$${\tilde a}_2\ (\ k_{\alpha}^{\prime},\ k_{\theta}^{\prime}\ )... ... 1.\ k_{\alpha}^{\prime} - 0.669716\ k_{\theta}^{\prime}\ )\ ;$$

$${\tilde a}_3\ (\ k_{\alpha}^{\prime},\ k_{\theta}^{\prime}\ )... ...6\ k_{\theta}^{\prime}\ )\ ; \ \ \ {\tilde a}_4\ =\ 16809.7\ ;$$

$$\Delta\ (\ k_{\alpha}^{\prime},\ k_{\theta}^{\prime}\ )\ = \ ... ...1952\ {k_{\alpha}^{\prime}}^2 - 1.\ {k_{\alpha}^{\prime}}^3\ -$$

$$-\ 0.0349361\ k_{\theta}^{\prime}\ +\ 0.541983\ k_{\alpha}^{\... ... \ -\ 2.03313\ {k_{\alpha}^{\prime}}^2\ k_{\theta}^{\prime}\ +$$

$$+\ 0.182602\ {k_{\theta}^{\prime}}^2\ -\ 1.37769\ k_{\alpha}^... ..._{\theta}^{\prime}}^2\ -\ 0.311139\ {k_{\theta}^{\prime}}^3\ )$$

and the output of the program

$$Region[\ \{\ \Delta > 0,\ a_0 > 0,\ a_1 > 0,\ a_2 > 0,\ a_3 >... ...rime},\ 0,\ 1.2\ \},\ \{\ k_{\alpha}^{\prime},\ 0,\ -1\ \}\ ]$$

will be the solution (domain) :

Figure 2.

On this picture it can be seen that asymptotical stability for given degree of stability could be reached by adding an accelerating dissipative force to the system, because $\ k_{\alpha}^{\prime} <0\ .$

REFERENCES

1. Banshchikov A.V., Bourlakova L.A. Information and Research System "Stability". //Informations RAS. The theory and control systems. (Izvestia RAN. Teoria i sistemi upravlenia.) $\ 1996.\ N^ {0} \ 2.\ 13 - 20.$ (in Russian).

2. Banshchikov A.V., Bourlakova L.A. On algorithms of symbolic computation at the research of the stability. // Programming. (Programirovanie.) $\ 1997.\ N^ {0} \ 3,\ 72-80.$ (in Russian).

3. Meierov M.V. Criterion of an aperiodicity of regulation. // Izvestia AN SSSR. Otdelenie tehnicheskix nauk. 1945. $\ N^{0}\ 12,\ 1169-1178.$ (in Russian). (Bulletin DE L'ACADÉMIE DES SCIENSES DE L'URSS. Classe des sciences techniques. )

4. Lankaster P. Theory of matrices. Academic Press. New York - London, 1969.

5. Potapenko E.M. Robust stability of controlled dynamical systems. // Izvestia RAN. Mehanika tverdogo tela. $\ 1993.\ V.5,\ 15-22.$ (in Russian).

6. Wolfram S. Mathematica: A System for Doing Mathematics by Computer. Addison-Wesley Publ.Co. 1988.

7. Potapenko E.M. Dynamics of a spacecraft with line active control by a gravitational stabilizer. // Space researches. (Kosmicheskie issledovania.) $\ 1988.\ V.\ 26.\ N^{0}\ 5,\ 699-708.$ (in Russian).