Optimization of Parameters of Asymptotically Stable Systems

Abstract

This work deals with numerical methods of parameter optimization for asymptotically stable systems. We formulate a special mathematical programming problem that allows us to determine optimal parameters of a stabilizer. This problem involves solutions to a differential equation. We show how to chose the mesh in order to obtain discrete problem guaranteeing the necessary accuracy. The developed methodology is illustrated by an example concerning optimization of parameters for a satellite stabilization system.

1. Introduction

Consider differential equation 𝑥̇=𝑓(𝑥,𝑢),𝑥∈𝑅^𝑛,𝑡≥0,(1.1) that describes a system equipped with a stabilizer. Here, 𝑢∈𝑈⊂𝑅^𝑘 is a parameter. It is assumed that 0=𝑓(0,𝑢) for all 𝑢∈𝑈 and the zero equilibrium position of system (1.1) is asymptotically stable whenever 𝑢∈𝑈. The parameter 𝑢 should be chosen to optimize, in some sense, the behaviour of the trajectories. The choice of this parameter can be based on various criteria; obviously, it is impossible to construct a stabilizer optimal in all aspects. For example, for a linear controllable system, the pole assignment theorem guarantees the existence of a linear feedback yielding a linear differential equation with any given set of eigenvalues. One can choose a stabilizer with a very high damping speed. However, such a stabilizer is practically useless because of the so called pick-effect (see [1, 2]). Namely, there exists a large deviation of the solutions from the equilibrium position at the beginning of the stabilization process, whenever the module of the eigenvalues is big.

The aim of this paper is to develop an effective numerical tool oriented to optimization of stabilizer parameters according to different criteria that appear in the engineering practice.

Throughout this paper, we denote the set of real numbers by 𝑅 and the usual 𝑛-dimensional space of vectors with components in 𝑅 by 𝑅^𝑛. The space of absolutely continuous functions defined in [0,𝑇] with values in 𝑅^𝑛 is denoted by AC([0,𝑇],𝑅^𝑛). We denote by ⟨𝑎,𝑏⟩ the usual scalar product in 𝑅^𝑛 and by |⋅| the Euclidean norm. By 𝐵, we denote the closed unit ball, that is, the set of vectors 𝑥∈𝑅^𝑛 satisfying |𝑥|≤1. The transpose of a matrix 𝐴 is denoted by 𝐴^∗. We use the matrix norm |𝐴|=max_|𝑥|=1|𝐴𝑥|. If 𝑃 and 𝑄 are two subsets in 𝑅^𝑛 and 𝜆∈𝑅, we use the following notations: 𝜆𝑃={𝜆𝑝∣𝑝∈𝑃}, 𝑃+𝑄={𝑝+𝑞∣𝑝∈𝑃,𝑞∈𝑄}.

2. Statement of the Problem

Denote by 𝑥(𝑡,𝑥₀,𝑢) the solution to the Cauchy problem 𝑥̇=𝑓(𝑥,𝑢),𝑥(0𝑥)∈=𝑅𝑥^𝑛₀,𝑡∈[0,𝑇],,(2.1) where 𝑢 is a parameter from a compact set 𝑈⊂𝑅^𝑘. Let 𝑓(0,𝑢)=0 for all 𝑢∈𝑈. Consider the functions 𝜑_𝑖(𝑢)=max_{𝑡∈Δ𝑖}max_{𝑥0∈𝐵𝑖}||𝑥𝑡,𝑥₀,𝑢||_𝑖,𝑖=0,𝑚.(2.2) Here, Δ_𝑖⊆[0,𝑇] are compact sets, and |⋅|_𝑖 are norms in 𝑅^𝑛, and 𝐵_𝑖={𝑥∈𝑅^𝑛∣|𝑥|_𝑖≤𝑏_𝑖}. Consider the following mathematical programming problem: 𝜑𝜑₀(𝑢)⟶min,_𝑖(𝑢)≤𝑢𝜑∈_𝑖,𝑈.𝑖=1,𝑚,(2.3) Many problems of stabilization systems' parameters optimization can be written in this form.

Minimization of the Final Deviation
The problem is to determine the optimal values of the system parameters that guarantee minimal deviation of the system state from the zero equilibrium position at the final moment of time. This problem can be formalized as follows: max_𝑥0∈𝐵||𝑥𝑇,𝑢𝑥₀∈,𝑢𝑈.||⟶min,(2.4) For linear systems 𝑥̇=𝐴(𝑢)𝑥 with 𝑇≫1, this problem is an approximation for the maximization of the degree of stability [3].

Minimization of the Maximal Deviation
This problem consists of determination of parameters that correspond to minimization of the maximum deviation of trajectories and satisfy certain restrictions at the final moment of time. This problem can be formalized as follows: max_{[]𝑡∈0,𝑇}mmaax_||𝑥0||=1x_||𝑥0||=1||𝑥𝑡,𝑥₀,𝑢||⟶min,||𝑥𝑢𝑇∈,𝑥₀𝑈,.𝑢||≤𝛿,(2.5)
The above problems are of interest for stabilization theory; they both have form (2.3). Problem (2.3) has some special features, and its solution can be useful for parameter optimization of stabilization systems; however, its study can hardly be performed analytically for more or less complex systems. For this reason, we focus on the numerical aspects of this problem.

3. Numerical Methods

Let 𝜀>0 be small enough. We approximate problem (2.3) by the following problem |||𝑥̃𝑡_𝑖𝑘𝜑₀⟶min,,𝑥^𝑖𝑗,𝑢|||_𝑖𝑢≤∈𝜑𝑈_𝑖+,𝜀,𝑖=0,𝑚,(3.1) where 𝑡^𝑖0=0, 𝑡^𝑖𝑘∈Δ_𝑖, 𝑥^𝑖𝑗∈𝐵_𝑖, 𝑗=1,𝐽, and 𝑥̃𝑡^𝑖𝑘+1,𝑥^𝑖𝑗,𝑢=𝑥̃𝑡^𝑖𝑘,𝑥^𝑖𝑗,𝑢+𝜏𝑓𝑥̃𝑡^𝑖𝑘,𝑥^𝑖𝑗,𝑢,𝑢,𝜏=𝑡^𝑖𝑘+1−𝑡^𝑖𝑘,𝑘=0,𝑁(3.2) is the Euler approximation for the solution 𝑥(⋅,𝑥^𝑖𝑗,𝑢). Problem (2.3) can be approximated by problems (3.1) with any given accuracy.

Assume that 𝑓(𝑥,𝑢)=𝐴(𝑢)𝑥+𝑔(𝑥,𝑢),(3.3) where matrix 𝐴(𝑢)=∇_𝑥𝑓(0,𝑢) has eigenvalues with negative real part and the function 𝑔(⋅,𝑢) satisfies 𝑔(0,𝑢)=0 and the Lipschitz condition ||𝑔𝑥₁,𝑢−𝑔𝑥₂,𝑢||≤𝐿^𝑢𝑔max||𝑥₁||,||𝑥₂||||𝑥₁−𝑥₂||,(3.4) with 𝐿^𝑢𝑔>0 for all 𝑥₁ and 𝑥₂ in a neighbourhood of the zero equilibrium position. Consider functions 𝜑_𝑖(⋅) defined by (2.2), assuming that the balls 𝐵_𝑖 are contained in a sufficiently small neighbourhood of the origin. Consider 𝛿>0. Let 𝐾_𝑖(𝛿) and 𝐽_𝑖(𝛿) be sets of indices such that the points 𝑡^𝑖𝑘∈Δ_𝑖, 𝑘∈𝐾_𝑖(𝛿), and 𝑥^𝑖𝑗∈𝐵_𝑖, 𝑗∈𝐽_𝑖(𝛿) form a 𝛿-net in Δ_𝑖 and 𝐵_𝑖, 𝑖=1,𝑚, respectively. Define the functions 𝜑^𝛿𝑖(𝑢)=max_{𝑘∈𝐾𝑖(𝛿)}max_{𝑗∈𝐽𝑖(𝛿)}|||𝑥̃𝑡_𝑖𝑘,𝑥^𝑖𝑗,𝑢|||,𝑖=0,𝑚.(3.5) Problem (3.1) can be written as 𝜑𝜑^𝛿0(𝑢)⟶min,_𝛿𝑖(𝑢)≤𝜑_𝑖+𝜀,𝑖=1,𝑚,𝑢∈𝑈.(3.6) Denote by 𝑢̂ and 𝑢^𝛿 the optimal parameters for problems (2.3) and (3.6), respectively.

Theorem 3.1. For any 𝜖>0, there exists 𝛿>0 such that 𝑢^𝛿 is an admissible solution to the following problem: 𝜑𝜑₀(𝑢)⟶min,_𝑖(𝑢)≤𝜑𝜑_𝑖+2𝜀,𝑖=1,𝑚,𝑢∈𝑈,₀𝑢^𝛿≤𝜑₀(𝑢)̂+2𝜀.(3.7)

This theorem allows one to choose the parameters of discretization in order to obtain optimal stabilizer parameters with a necessary precision. A rigorous formulation of this claim is the following. Denote by 𝑉(𝜎) the value of the problem 𝜑𝜑₀(𝑢)⟶min,_𝑖(𝑢)≤𝜑_𝑖+𝜎,𝑖=1,𝑚,𝑢∈𝑈.(3.8) Assume that problem (2.3) is calm in Clarke's sense (see [4]). Then, there exists a constant 𝐾>0 satisfying the inequality 𝑉(2𝜀)−𝑉(0)2𝜀>−𝐾,(3.9) for all 𝜀>0 sufficiently small. It follows from Theorem 3.1 that ||𝑉−𝑉(0)||≤𝑀𝜀,(3.10) where 𝑉=𝜑₀(𝑢)̃, 𝑢̃ is the solution of problem (3.1), and 𝑀=2max{1,𝐾}.

The exact formulas for 𝛿=𝛿(𝜖) leading to the proof of Theorem 3.1 are presented in the Appendix. The main tool used to obtain them is the following version of Filippov-Gronwall inequality [5].

Theorem 3.2. Let 𝑃={𝑝∈𝑅^𝑛∣⟨𝑝,𝑉𝑝⟩≤1}, where 𝑉 is a symmetric positive definite matrix. Consider the functions 𝑦(⋅)∈AC([0,𝑇],𝑅^𝑛) and 𝜉(⋅)∈AC([0,𝑇],𝑅), 𝜉(𝑡)≥0 satisfying the following condition max_{⟨𝑝,𝑉𝑝⟩=1}(⟨𝑦̇(𝑡),𝑉𝑝⟩−⟨𝑓(𝑦(𝑡)−𝜉(𝑡)𝑝),𝑉𝑝⟩)≤𝜉̇(𝑡),(3.11) for almost all 𝑡∈[0,𝑇]. Then, 𝑥(𝑡)∈𝑦(𝑡)+𝜉(𝑡)𝑃 for all 𝑡∈[0,𝑇], whenever 𝑥₀∈𝑦(0)+𝜉(0)𝑃, where 𝑥(𝑡) is the solution to the Cauchy problem 𝑥̇=𝑓(𝑥), 𝑥(0)=𝑥₀.

Note that the use of this theorem allows us to obtain more precise estimates for the number of points in the meshes needed to achieve a given discretization accuracy. The estimates based on the usual Gronwall inequality can be significantly improved for asymptotically stable systems if we take into account the behaviour of the trajectories for large values of time. Theorem 3.2 is a natural tool for this analysis. For example, according to the classical estimates, the number of points in the mesh in 𝑡, needed to ensure a given precision, grows exponentially with the length of the time interval. Meanwhile, the estimates obtained from Theorem 3.2 for asymptotically stable systems (see Theorems A.2 and A.6) give a linear growth of the number of points in the mesh. This result is of practical importance. Optimization problem (3.6) is a hard nonsmooth problem. Our computational experience shows that the usage of the Nelder-Mead method is the most adequate approach to solve it. The numerical solution of this problem significantly depends on the structure of the involved functions. The problem of optimal choice of parameters is solved only once, at the stage of the control system's development, so one could afford to dedicate more resources to its solution. However, if the mesh is constructed using the classical precision estimates, the required computational effort can be extremely high, making it impossible to solve the problem in a reasonable time. Our estimates for the number of points of discretization allow us to construct an adequate mesh and to significantly reduce the CPU time.

4. Example: Optimal Parameters for Satellite-Stabilizer System

Consider motion of a connected two-body system in a circular orbit around the Earth. Body 1 is a satellite with the center of mass 𝑂₁, and body 2 is a stabilizer with the center of mass 𝑂₂. These two bodies are linked to each other at the point 𝑃 through a dissipative hinge mechanism (Figure 1). Let 𝑂 be the center of mass of the system.

Figure 1 

The satellite stabilizer system.

We use three reference frames: 𝑂𝑋𝑌𝑍 is the orbital coordinate frame, its axis 𝑂𝑍 is directed along the radius vector of the point 𝑂 with respect to the center of the Earth, 𝑂𝑋 is directed along the velocity of the point 𝑂, and 𝑂𝑌 is normal to the orbit plane. The axes of referential frames 𝑂₁𝑥₁𝑦₁𝑧₁ and 𝑂₂𝑥₂𝑦₂𝑧₂ are the central principal axes of inertia for bodies 1 and 2, respectively. Consider motion of the system in the orbit plane supposing that the bodies are connected in their centres of mass; that is, the points 𝑂₁, 𝑂₂, 𝑂, and 𝑃 coincide. Let 𝛼₁ and 𝛼₂ be the angles between the axis 𝑂𝑋 and the axes 𝑂₁𝑥₁ and 𝑂₂𝑥₂, respectively. Denote by 𝛼^𝑖, 𝑖=1,2 the derivative of 𝛼_𝑖 with respect to time. The equations of motion for this system can be written as [6]𝐵𝐵₁𝛼^1+3𝜔²⁰𝐴₁−𝐶₁sin𝛼₁cos𝛼₁+𝑘₁𝛼^1−𝛼^2=0,₂𝛼^2+3𝜔²⁰𝐴₂−𝐶₂sin𝛼₂cos𝛼₂−𝑘₁𝛼^1−𝛼^2=0.(4.1) Here, 𝐴₁, 𝐵₁, 𝐶₁ and 𝐴₂, 𝐵₂, 𝐶₂ are the principal moments of inertia of the bodies,𝑘₁ is the damping coefficient of the system, and 𝜔₀ is the constant angular velocity of the orbital motion of the system's center of mass. Introducing a new independent variable 𝜏=𝜔₀𝑡 and the dimensionless parameters 𝑝₁=𝐴¹−𝐶₁𝐵₁,𝑝₂=𝐴²−𝐶₂𝐵₂,𝜇=𝐵²𝐵₁,𝑘₁=𝑘₁𝜔₀𝐵₁,(4.2) the equations of motion can be written as 𝛼𝛼̈₂̈₁+3𝑝₁sin𝛼₁cos𝛼₁+𝑘₁𝛼̇₁−𝛼̇₂=0,+3𝑝₂sin𝛼₂cos𝛼₂−𝑘¹𝜇𝛼̇¹−𝛼̇₂=0.(4.3) Here, the dot denotes the derivative with respect to 𝜏. The parameters (𝑝₁,𝑝₂,𝑘₁,𝜇) satisfy the following conditions: −1≤𝑝₁≤1,−1≤𝑝₂≤1,𝜇>0,𝑘₁>0.(4.4) We study small oscillations of system (4.3) in the vicinity of the equilibrium position 𝛼₁₀=0,𝛼₂₀=0.(4.5) The equations of motion, linearized in the vicinity of the above stationary solution, take the form 𝛼𝛼̈₂̈₁+3𝑝₁𝛼₁+𝑘₁𝛼̇₁−𝛼̇₂=0,+3𝑝₂𝛼₂−𝑘¹𝜇𝛼̇¹−𝛼̇₂=0.(4.6) The characteristic equation for system (4.6) is 𝜇𝜆⁴+𝑘₁(1+𝜇)𝜆³+3𝜇𝑝₁+𝑝₂𝜆²+3𝑘₁𝑝₁+𝜇𝑝₂𝜆+9𝜇𝑝₁𝑝₂=0.(4.7) Analysis of (4.7) allows one to obtain the necessary and sufficient conditions of asymptotic stability. The region of asymptotic stability is given by 𝑘₁,𝑝₁,𝑝₂∶𝑘₁>0,𝑝₁>0,𝑝₂>0,𝑝₁≠𝑝₂.(4.8) Taking into account the feasibility conditions for the system parameters, we arrive at the following set of admissible parameters for our optimization problem: 𝑈=𝑝₁,𝑝₂,𝑘₁,𝜇∶𝑘₁>0,𝜇>0,0<𝑝₁≤1,0<𝑝₂≤1,𝑝₁≠𝑝₂.(4.9)

4.1. The Maximal Degree of Stability

Consider the set 𝑈 described by (4.9). Denote by 𝑢=(𝑝₁,𝑝₂,𝑘₁,𝜇) a parameter that belongs to the set 𝑈. Let {𝜆₁(𝑢),…,𝜆₄(𝑢)} be the roots of (4.7). The inclusion 𝑢∈𝑈 implies that Re𝜆_𝑖(𝑢)<0, 𝑖=1,4. The degree of stability is defined by 𝛿(𝑢)=−max_𝑖=1,4Re𝜆_𝑖(𝑢).(4.10) Consider the following problem 𝛿(𝑢𝑢)⟶∈𝑈m.ax,(4.11) In [7, 8], it is proved that the maximal degree of stability is achieved when all the roots of the characteristic equations are real and equal. This situation becomes possible only when the conditions 𝑘₁=4𝛿(𝑢)𝜇𝑝1+𝜇₁,3𝑝+𝑝₂=2𝛿²(𝑢),₁+𝜇9𝑝𝑝₂₁=(1+𝜇)𝛿²(𝑢),𝑝₂=𝛿⁴(𝑢)(4.12) are satisfied. The above system has two sets of solutions 𝑝₁̂=3−𝑝2₂√̂2²≃0.0294,𝑘̂=1,₁=𝜇√=63−32−√2√2≃0.4203,𝑝₁̂2≃0.1716,(4.13)𝑝₂̂=1,=𝑘̂3−2√2²≃0.0294,₁𝜇=3=+√2√6≃2.4495,2≃5.8284.(4.14)

4.2. Numerical Optimization

Denote by 𝑥(⋅,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇) the solution of linear system (4.6) with 𝑥(0)=𝑥₀, defined in the interval [0,𝑇]. The parameters (𝑝₁,𝑝₂,𝑘₁,𝜇) belong to asymptotic stability region defined in (4.9). Consider the following problem: max_||𝑥0||=1||𝑥𝑇,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇||⟶min,𝑘₁>0,𝜇>0,0<𝑝₁≤1,0<𝑝₂≤1.(4.15) Problem (4.15) can be reduced to an optimization problem without constrains using quadratic penalty functions; see [9]. If 𝑇 is big enough, this problem approximates problem (4.11), where the concept of degree of stability is used. The parameters given by (4.13) and (4.14) are close to optimal solutions of problem (4.15) only when 𝑇≫1. Put 𝑇=10𝜋. The results of simulations show that the value of problem (4.15) is about 10⁻⁶-10⁻⁷, independently on the values of admissible parameters (𝑝₁,𝑝₂,𝑘₁,𝜇). For example, the following values 𝑝̂₁̂=0.0779,𝑝̂₂̂=0.8574,𝑘̂̂₁=0.5540,𝜇₁=0.3337(4.16) are optimal parameters for problem (4.15). The corresponding minimal value is 7.2×10⁻⁷. Parameters (4.13) and (4.14) give the values 1.0×10⁻⁶ and 5.9×10⁻⁷, respectively.

In practice, it is important to consider smaller time intervals [0,𝑇]. Solve problem (4.15) with 𝑇=3𝜋. In this case, we see that the value of problem really depends on the choice of parameters. Let us estimate the global minimum in this problem. To this end, we solve problem (4.15) using all combinations of the following values: 𝑝𝑝₁=0.25,0.5,0.75,₂=𝑘0.25,0.5,0.75,₁𝜇==21,,42,,63,,(4.17) as initial guesses for numerical optimization. We obtain the following two sets of parameters with the best value of the problem: 𝑝𝑝₁₁=0.06928,𝑝₂=1.00757,𝑘₁=0.59209,𝜇=0.33161,(4.18)=1.00521,𝑝₂=0.06920,𝑘₁=1.78178,𝜇=3.01152.(4.19) The estimate for the global minimum 𝑚 is 𝑚=0.00378.(4.20) Meanwhile, the value corresponding to parameters (4.13) and (4.14) is 0.4660. Thus, we see that the methodology based on resolution of problem (2.3) can be more adequate in the practice than that one using the concept of degree of stability.

Since we are studying the behaviour of a nonlinear system in a vicinity of its equilibrium position, it is also important to estimate the deviation of the linearized system trajectories from zero. The stabilizer constructed for the linearized system makes sense only if its trajectories belong to a small vicinity of the equilibrium position; otherwise, the influence of the nonlinearity can destabilize the system even in a very small neighbourhood of the equilibrium.

Consider the following problem: mamx_{[]𝑡∈0,3𝜋}axm_||𝑥0||=1ax_||𝑥0||=1||𝑥𝑡,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇||⟶min,||𝑥3𝜋,𝑥₀,𝑘𝑝₁,𝑝₂,𝑘₁,𝜇||≤0.005,₁0<𝜇𝑝>>₁00,,0<𝑝₂≤1,≤1.(4.21) In this problem, the solutions of system (4.6) at the moment 𝑇=3𝜋 are constrained to be in a neighbourhood of the equilibrium position with the radius 0.005. The obtained optimal solutions (𝑝₁̃,𝑝₂̃,𝑘̃₁,𝜇) minimize the maximum norm of the damping process of linear system (4.6) in the interval [0,3𝜋]. After numerical optimization, we get the following optimal parameters: 𝑝𝑝₁₁̃̃=0.07140,𝑝₂̃=1.01643,𝑘̃₁=0.60004,𝜇=0.33887,(4.22)=1.01642,𝑝₂̃=0.07140,𝑘̃₁=1.77067,𝜇=2.95097.(4.23) The corresponding value of problem (4.21) is 𝑃=1.58685. We can see that the couple of parameters in (4.22) and (4.23) are slightly different from (4.18) and (4.19). Moreover,max_{[]𝑡∈0,3𝜋}max_||𝑥0||=1|||𝑥𝑡,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇|||≃1.5974≃𝑃.(4.24) Taking the optimal parameters (𝑝₁̂,𝑝₂̂,𝑘̂₁,𝜇) of problem (4.11), we get max_{𝑡∈[0,3𝜋]}max_||𝑥0||=1||𝑥𝑡,𝑥₀,𝑝₁̂,𝑝₂̂,𝑘̂₁,𝜇||≃1.9169.(4.25) Thus, the stabilizer with the parameters corresponding to the maximal degree of stability yields more significant deviation of the trajectories from the equilibrium position than the stabilizer with the parameters obtained solving problem (4.21).

Our aim is to find optimal parameters for system (4.3). To this end, we solve problem (4.15), with 𝑇=3𝜋, but now, 𝑥(⋅,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇) stands for the solution of system (4.3) with 𝑥(0)=𝑥₀. We get the following two sets of optimal parameters: 𝑝₁=0.23350,𝑝₂=1.08235,𝑘₁=0.62791,𝜇=0.62137,(4.26)𝑝₁=1.07743,𝑝₂=0.23171,𝑘₁=1.02413,𝜇=1.63140.(4.27) Consider optimal parameters (4.26) and (4.27), and denote by 𝑥(⋅,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇) the solution of system (4.3) corresponding to these parameters and satisfying 𝑥(0)=𝑥₀. The optimal value is 𝑃₁=max_{𝑡∈[0,3𝜋]}max_||𝑥0||=1|||𝑥𝑡,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇|||≃1.3634.(4.28) The above value of 𝑃₁ is even smaller than 𝑃, obtained for the linearized system.

To compare the solutions of all optimization problem, we consider the following function 𝑀(𝛿)=max_{||𝑥0||≤𝛿}||𝑥3𝜋,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇||,𝛿>0.(4.29) Here, 𝑥(3𝜋,𝑥₀,𝑝₁,𝑝₂,𝑘₁,𝜇) stands for the solution to system (4.3). Table 1 gives the values of the function 𝑀(𝛿) for parameters (4.13), (4.18), (4.22), and (4.26).

Table 1 

Values of function 𝑀(𝛿) corresponding to parameters (4.13), (4.18), (4.22), and (4.26).

Observe that the values in the last column of Table 1, computed for parameters (4.26), always satisfy the condition 𝑀(𝛿)≤0.1𝛿. On the other hand, this condition is not satisfied for parameters (4.13), obtained maximizing the degree of stability. For the parameters obtained for linearization, the condition is satisfied only if 𝛿 is sufficiently small. This illustrates the advantages of applying the introduced methodology, based on numerical solution of optimization problem (2.3), directly to nonlinear systems.

5. Conclusions

The methods usually applied to optimize the parameters of a stabilization system are based on the idea of the maximum stability degree, that is, the minimization of the system's transition time. These methods, however, face the problem of the so-called peak effect when the deviation of the system trajectory from the equilibrium increases with the decrease of the time of response. The approach suggested in this paper consists of a numerical analysis of a stabilization system based on a more complete and flexible description of the system behaviour capable to take into account not only the stability degree but also the maximum deviation of the trajectory on a given time interval or at a given moment. For this optimization problem, we develop a numerical method and prove that it can guarantee a given accuracy for the problem solution. We obtained more precise estimates for the number of points in the meshes needed to achieve a given discretization accuracy than the estimates based on the usual Gronwall inequality. This method is applied to optimization of a stabilization system for a satellite with a gravitational stabilizer. The obtained results show that the above approach can offer solutions more adequate for practical implementation than those given by optimization of the stability degree.

Appendix

A. The Mathematical Basis

In this Appendix, we present a series of theorems, with schematic proofs, containing explicit estimates for the fineness of discretization needed to obtain the necessary precision of approximations and to prove Theorem 3.1.

A.1. Linear Systems

Consider a linear system 𝑥̇=𝐴𝑥,𝑥∈𝑅^𝑛,𝑡≥0,(A.1) where 𝐴 is a matrix. Assume that all its eigenvalues have negative real part. Let 𝑉 be a symmetric positive definite matrix satisfying the Lyapunov equation [10] 𝐴^∗𝑉+𝑉𝐴=−𝐼.(A.2) Set 𝑃={𝑝∈𝑅^𝑛∣⟨𝑝,𝑉𝑝⟩≤1}.(A.3) The quadratic form 𝑉(𝑝)=⟨𝑝,𝑉𝑝⟩ is the Lyapunov function for system (A.1). Denote by 𝜂₁ and 𝜂₂ the minimal and maximal eigenvalues of 𝑉, respectively. Let 𝜏 be a positive constant. Consider the Euler approximation for system (A.1), 𝑦_𝑘+1=𝑦_𝑘+𝜏𝐴𝑦_𝑘,𝑘=0,1,….(A.4) The following theorem provides an explicit estimate for the constant 𝜏 guaranteeing the equality lim_𝑘→∞𝑦_𝑘=0.

Theorem A.1. Let 𝑏>0. Consider 𝑦₀∈𝑏𝑃. If 0<𝜏<𝜂¹𝜂²²||𝐴||²,(A.5) then the following inequalities hold: ⟨𝑦_𝑘,𝑉𝑦_𝑘⟩≤𝛽^𝑘⟨𝑦₀,𝑉𝑦₀⟩,𝑘=1,2,…,(A.6) where 𝛽=1−𝜏𝜂₂+𝜏²𝜂₂𝜂₁||𝐴||².(A.7)

It is easy to see that 0<𝛽<1. The proof of this theorem uses the induction and the inequality 𝜂₁||𝑝||²≤⟨𝑝,𝑉𝑝⟩≤𝜂₂||𝑝||².(A.8)

Assume that constant 𝜏 satisfies condition (A.5). Consider the polygonal Euler approximation to solution of system (A.1) 𝑦(𝑡)=𝑦_𝑘𝑡_𝑘+𝑡−𝑡_𝑘𝐴𝑦_𝑘,𝑡∈𝑡_𝑘,𝑡_𝑘+1,=𝑘𝜏,𝑘=0,1,….(A.9) Let 𝑏>0. Set 𝜉𝛾=−min_{⟨𝑝,𝑉𝑝⟩=1⟨𝑦,𝑉𝑦⟩≤𝑏2}⟨𝐴²𝑦,𝑉𝑝⟩,(A.10)_𝜇(𝑡)𝑥=(𝑡)𝜉=₀𝑒𝑒^−𝜇𝑡_𝐴𝑡,𝜉₀>0,(A.11)𝑥(0).(A.12)

Theorem A.2. Let 𝑦₀∈𝑏𝑃 be given. Assume that condition 0≤𝜇<1/2𝜂₂ is satisfied. If 𝜏𝛾≤12𝜂₂−𝜇𝜉_𝜇(𝑡),𝑡≥0,(A.13) then the inequality |𝑦(𝑡)−𝑥(𝑡)|≤𝜉_𝜇(𝑡)/√𝜂₁ holds for all 𝑡∈[0,𝑇], whenever |𝑦₀−𝑥(0)|≤𝜉₀/√𝜂₂.

This theorem is a consequence of the inequality max_{⟨𝑝,𝑉𝑝⟩=1}⟨𝐴𝑝,𝑉𝑝⟩≤−12𝜂₂,(A.14) of the inclusion 1√𝜂₂𝐵⊂𝑃⊂1√𝜂₁𝐵,(A.15) and of Theorem 3.2 with function 𝑦(⋅) defined by (A.9) and function 𝜉(⋅) defined by (A.11).

The following theorem is also a consequence of Theorem 3.2.

Theorem A.3. Let 𝑥(𝑡)=𝑒^𝐴𝑡𝑥₀ and 𝑧(𝑡)=𝑒^𝐴𝑡𝑧₀ be solutions to system (A.1). Assume that 𝑥₀,𝑧₀∈𝑃. If 0≤𝜇≤12𝜂₂,(A.16) then the inequality |𝑧(𝑡)−𝑥(𝑡)|≤𝜉_𝜇(𝑡)/√𝜂₁, 𝑡∈[0,𝑇] holds whenever |𝑧₀−𝑥₀|≤𝜉₀/√𝜂₂.

Consider now 𝑡∈Δ⊂[0,𝑇], where Δ is a closed interval. Let 𝑥(𝑡,𝑥₀) be the solution of system (A.1), with 𝑥₀∈𝑆={𝑥∈𝑅^𝑛∶|𝑥|=1}.(A.17) Let 𝛿>0. Assume that parameters of function 𝜉_𝜇(⋅) defined by (A.11) satisfy the following conditions, 𝜉₀=√𝜂₂𝛿,0≤𝜇<12𝜂₂.(A.18) Assume that {𝑥_𝑗}, 𝑗=1,𝐽 is a finite set of points uniformly distributed on 𝑆. If 𝐽≥21𝛿+1^𝑛−1,(A.19) then we have ^𝐽_𝑗=1𝑥_𝑗+𝛿𝐵⊃𝑆.(A.20) Let |Δ| be the length of the interval Δ. Consider a finite set {𝑡_𝑘}⊂Δ, 𝑘=0,𝑁, such that the difference 𝑡_𝑘+1−𝑡_𝑘=𝜏=𝛿2𝛾√𝜂₂(A.21) is a constant. If 𝑁=2𝛾√𝜂₂𝛿||Δ||+1,(A.22) then the set ^𝑁_𝑘=0𝑡^𝑘−𝛿4𝛾√𝜂₂,𝑡_𝑘+𝛿4𝛾√𝜂₂(A.23) contains Δ. Consider the Euler approximation of solution 𝑥(𝑡,𝑥₀)𝑥̃𝑡_𝑘+1,𝑥_𝑗=𝑥̃𝑡_𝑘,𝑥_𝑗+𝜏𝐴𝑥̃𝑡_𝑘,𝑥_𝑗,𝑘=0,𝑁.(A.24)

Theorem A.4. Let 𝜀>0. If 𝛿=4𝛾√𝜂₁8𝛾√𝜂₂+||𝐴||𝜀,(A.25) then the following inequality holds: ||||max_𝑘=0,𝑁max_𝑗=1,𝐽||𝑥̃𝑡_𝑘,𝑥_𝑗||−max_𝑡∈Δmax_𝑥0∈𝑆||𝑥𝑡,𝑥₀||||||≤𝜀.(A.26)

The proof of this theorem uses the results of Theorems A.1, A.2, and A.3.

A.2. Nonlinear Systems

Assume that 𝑔∶𝑅^𝑛→𝑅^𝑛 is a twice continuously differentiable function satisfying 𝑔(0)=0. Consider the function 𝑓(𝑥)=𝐴𝑥+𝑔(𝑥), where 𝐴 is a matrix. Assume that the eigenvalues of the matrix 𝐴 have negative real parts. Consider the system 𝑥̇=𝑓(𝑥),𝑡≥0.(A.27) Since 𝑔 is twice continuously differentiable, there exists a constant 𝐿_𝑔>0 such that function 𝑔(⋅) satisfies the following Lipschitz condition: ||𝑔𝑥₁−𝑔𝑥₂||≤𝐿_𝑔max||𝑥₁||,||𝑥₂||||𝑥₁−𝑥₂||,(A.28) for all 𝑥₁ and 𝑥₂ in a small neighbourhood of the equilibrium position 𝑥=0. Consider the set 𝑃 defined by (A.3) and constants 𝜂₁, 𝜂₂ as before. Define the Euler approximation for system (A.27), 𝑦_𝑘+1=𝑦_𝑘+𝜏𝑓𝑦_𝑘,𝑘=0,1,…,(A.29) where 𝜏 is a positive constant.

Theorem A.5. Assume that 𝑏 is a constant satisfying 0<𝑏<√𝜂₁4𝐿_𝑔𝜂₂.(A.30) Let 𝑦₀∈𝑏𝑃. If 0<𝜏<8𝜂²¹𝜂₂4||𝐴||𝜂₂+1²,(A.31) then the following inequalities hold: ⟨𝑦_𝑘,𝑉𝑦_𝑘⟩≤𝛽^𝑘⟨𝑦₀,𝑉𝑦₀⟩,𝑘=1,2…,(A.32) where 𝛽=1−𝜏2𝜂₂+𝜏²4||𝐴||𝜂₂+1²16𝜂²¹.(A.33)

The proof of this theorem uses the Lipschitz condition (A.28) in the form ||𝑔𝑦_𝑘||≤𝐿_𝑔||𝑦_𝑘||²,(A.34) and the mathematical induction method. It is easy to see that 0<𝛽<1.

Assume that the constant 𝜏 satisfies condition (A.31) and consider the function 𝑦(𝑡)=𝑦_𝑘+𝑡−𝑡_𝑘𝑓𝑦_𝑘,𝑡∈𝑡_𝑘,𝑡_𝑘+1,𝑡_𝑘=𝑘𝜏,𝑘=0,1,….(A.35) Put 𝛾₁=𝛾+8𝜂²||𝐴||+164𝜂²²√𝜂₁𝐿_𝑔,(A.36) where 𝛾 is defined by (A.10). Denote by 𝑥(𝑡) the solution to system (A.27) with the initial condition 𝑥(0)=𝑥₀.

Theorem A.6. Assume that 𝑏 is a constant satisfying 𝜉₀≤𝑏≤√𝜂₁4𝐿_𝑔𝜂₂,(A.37) and let 𝑦₀∈𝑏𝑃. If 0≤𝜇<1/4𝜂₂ and 𝜏𝛾₁≤14𝜂₂−𝜇𝜉_𝜇(𝑡),𝑡≥0,(A.38) then the inequality |𝑦(𝑡)−𝑥(𝑡)|≤𝜉_𝜇(𝑡)/√𝜂₁ holds for all 𝑡∈[0,𝑇], whenever |𝑦₀−𝑥(0)|≤𝜉₀/√𝜂₂.

This theorem is a consequence of Theorem 3.2 with function 𝑦(⋅) defined by (A.35) and function 𝜉(⋅) defined by (A.11).

Consider now the Cauchy problems 𝑥𝑥̇(=0)𝑓=(𝑥𝑥)₀,𝑧𝑧̇(=0)𝑓=(𝑧𝑧)₀,,.(A.39)

Denote by 𝑥(⋅) and 𝑧(⋅) the respective solutions. The following theorem is also a consequence of Theorem 3.2.

Theorem A.7. Let 𝑏 be a constant satisfying 𝜉₀≤𝑏≤√𝜂₁4𝐿_𝑔𝜂₂.(A.40) Assume that the trajectories 𝑥(⋅) and 𝑧(⋅) belong to the set 𝑏𝑃. If 0≤𝜇≤14𝜂₂,(A.41) then we have |𝑧(𝑡)−𝑥(𝑡)|≤𝜉_𝜇(𝑡)/√𝜂₁ for all 𝑡∈[0,𝑇], whenever |𝑧₀−𝑥₀|≤𝜉₀/√𝜂₂.

Let 𝛿>0. Assume that the parameters of function 𝜉_𝜇(⋅) defined by (A.11) satisfy the following conditions: 𝜉₀=√𝜂₂𝛿,0≤𝜇≤14𝜂₂.(A.42) Consider the balls𝐵_𝑏𝑖=𝑥∈𝑅^𝑛∶|𝑥|≤𝑏^𝑖√𝜂₂⊂𝑏^𝑖𝑃,𝑖=0,𝑚,(A.43) where the constants 𝑏_𝑖 satisfy the following conditions: 𝜉₀≤𝑏_𝑖≤√𝜂₁4𝐿_𝑔𝜂₂,𝑖=0,𝑚.(A.44) For each index 𝑖, take a finite set {𝑥^𝑖𝑗} of points, 𝑗=1,𝐽_𝑖, uniformly distributed in the ball 𝐵_𝑏𝑖. If 𝐽_𝑖≥𝑏^𝑖√𝜂₂𝛿+1^𝑛,(A.45) then we have ^𝐽𝑖_𝑗=1𝑥^𝑖𝑗+𝛿𝐵⊃𝐵_𝑏𝑖.(A.46) Let Δ_𝑖⊂[0,𝑇], 𝑖=0,𝑚, be closed intervals with length |Δ_𝑖|. Consider a finite set of points {𝑡^𝑖𝑘}, 𝑘=0,𝑁_𝑖, in each interval Δ_𝑖. It is assumed that the difference 𝜏=𝑡^𝑖𝑘+1−𝑡^𝑖𝑘 is a constant. Let 𝜏=𝛿4𝛾₁√𝜂₂.(A.47) Define the sets Δ_𝑁𝑖=^𝑁𝑖_𝑘=0𝑡^𝑖𝑘−𝛿8𝛾₁√𝜂₂,𝑡^𝑖𝑘+𝛿8𝛾₁√𝜂₂,𝑖=0,𝑚.(A.48) If 𝑁_𝑖=4𝛾¹√𝜂₂𝛿||Δ^𝑖||+1,(A.49) then we have Δ_𝑁𝑖⊃Δ_𝑖. Let 𝑥₀∈𝐵_𝑏𝑖. Denote by 𝑥(𝑡,𝑥₀) the solution to the Cauchy problem 𝑥̇=𝐴𝑥+𝑔𝑥(𝑥()0,)=𝑥𝑥∈₀𝑅^𝑛,𝑡∈Δ_𝑖,.(A.50) Consider the Euler approximation of the solution 𝑥(𝑡,𝑥₀)𝑥̃𝑡^𝑖𝑘+1,𝑥^𝑖𝑗=𝑥̃𝑡^𝑖𝑘,𝑥^𝑖𝑗+𝜏𝐴𝑥̃𝑡^𝑖𝑘,𝑥^𝑖𝑗+𝑔𝑥̃𝑡^𝑖𝑘,𝑥^𝑖𝑗,𝑘=0,𝑁_𝑖,(A.51) with 𝜏 satisfying condition (A.47).

Theorem A.8. Let 𝜀>0 be given. If 𝛿=2⁷𝐿_𝑔𝛾₁𝜂²²√𝜂₁𝜂₂2⁸𝐿_𝑔𝛾₁𝜂³²+4𝜂₂√𝜂₁||𝐴||+√𝜂₁𝜀,(A.52) then the following inequalities: ||||max_{𝑘=0,𝑁𝑖}max_{𝑗=1,𝐽𝑖}|||𝑥̃𝑡_𝑖𝑘,𝑥^𝑖𝑗|||−max_{𝑡∈Δ𝑖}max_{𝑥0∈𝐵𝑏𝑖}||𝑥𝑡,𝑥₀||||||≤𝜀,𝑖=0,𝑚,(A.53) hold.

The proof of this theorem follows from Theorems A.5, A.6, and A.7.

Acknowledgments

The authors are grateful to the reviewers for their valuable suggestions. This research is supported by the Portuguese Foundation for Science and Technologies (FCT), the Portuguese Operational Programme for Competitiveness Factors (COMPETE), the Portuguese Strategic Reference Framework (QREN), and the European Regional Development Fund (FEDER).