Synergetic Control Theory Approach for Solving Systems of Nonlinear Equations Anton BEZUGLOV Earth Sciences Resources Institute, University of South Carolina Columbia, SC, 29208 Anatoliy KOLESNIKOV Synergetics and Controlling Processes Dept, Taganrog State University of RadioEngineering, Taganrog, Russian Federation Igor KONDRATIEV Electrical Engineering Dept, University of South Carolina, Columbia, SC, 29208 Juan VARGAS Computer Science & Engineering Dept, University of South Carolina, Columbia, SC, 29208
ABSTRACT This paper presents Synergetic Control Theory (SCT) and discusses how it can be used for solving systems of nonlinear equations. SCT is a new methodology for solving systems of nonlinear equations. The main advantage of SCT is that it maps the original system of equations to a dynamical system such that (1) any trajectory in the state space of the system ends in an attracting point; (2) the attracting point is located at the solution of the original system and (3) the rate at which the dynamical system moves towards the attracting point is controllable. An algorithm based on SCT is discussed. The algorithm has the following advantages: (1) If solutions exist, the algorithm finds one at a controllable rate, independently of initial guesses, (2) If no solutions exist, the algorithm makes such determination quickly, using stability analysis. Keywords: Nonlinear equations, Systems of nonlinear equations, dynamical systems, stability analysis, Synergetic Control Theory, attracting regions, attractors. 1. INTRODUCTION The purpose of this work is to describe a synergetic method for obtaining solutions to systems of nonlinear equations. Similar to the traditional Newton-Raphson and Broyden’s methods, the synergetic method transforms the equations to a dynamical system that performs iterations. The synergetic dynamical system can be considered superior to other systems, because: (1) It is linear and (2) It has attractors only where they are needed, i.e. at the solutions. This paper is organized as follows. Section 2 discusses the Newton-Raphson and Broyden’s methods. In section 3 the basics of SCT are given and the
synergetic algorithm is explained. Section 4 provides examples of finding roots using the synergetic algorithm and compares the convergence rate versus the NewtonRaphson method. 2. BACKGROUND Iterative methods for finding roots of nonlinear systems of equations are used, because mathematics does not have a sufficiently general analytical method. In contrast to analytical methods, iterative methods consider the process of solving equations or systems of equations as a dynamical process, which controls the convergence from an initial guess to a solution. The dynamical systems that guide these processes are a by product of the iterative methods. Therefore, the properties of the methods, such as convergence rate, robustness to initial guesses, robustness to local maxima and minima, are mostly determined by the properties of the dynamical systems. A typical approach is the expansion of the Taylor’s series, which is applied in a family of methods, including the Newton-Raphson, the Broyden’s, and their modifications [1,2]. Let us review approaches used in these methods. Suppose that a problem is formulated as follows:
Fi (x ) = 0
(1)
where Fi , 1 ≤ i ≤ M is a set of equations that need to be zeroed, x is a vector of variables. The vector of functions F can be approximated as a part of the Taylor’s series as:
F ( x + ∆x ) = F ( x ) + J ⋅ ∆x
(2)
where J is the Jacobian. The iterative formulae can be
derived by setting F ( x + ∆x ) = 0 and solving (2) for
∆x :
∆x = − J −1 ⋅ F ( x )
(3)
The Newton-Raphson method The Newton-Raphson method uses the dynamical system based on equation (3). The solution is approached from an initial guess x0 by the iterative formula:
xk +1 = xk + ∆x
(4)
where k is an iteration number. The converging process stops when the distance iterations become less than some threshold: xk +1 − xk < ε , ε > 0 .
Figure 1. Failures of the Newton-Raphson method
The Newton-Raphson method is known to have a quadratic convergence rate, provided that the initial guesses happen to be near the roots. If the initial guesses are not near the roots, the method may fail to converge. Figure 1 (taken from [3]) shows several types of unsuccessful choices of initial guesses. In the top left, the Newton-Raphson will infinitely loop at the local minimum. The top right and bottom right demonstrate a numerical overflow if the guess is close to a point where the first derivative of F (⋅) is close to zero. In the bottom left graph, the divergence is demonstrated. The insufficient global convergence and instability at the local extrema are limitations of the NewtonRaphson method, which can be improved. For instance, the iteration formula can be modified xk +1 = xk + α ⋅ ∆x , where 0 < α ≤ 1 . The parameter is chosen so that it minimizes F ( xk +1 ) . This modification, which is supposed to enhance the convergence, is a heuristic, and as such, does not guarantee a global convergence.
α
2
The Broyden’s method The Broyden’s method [4] does not require an exact Jacobian, this may save significant computational time as compared to the Newton-Raphson method. The Broyden’s method is a secant method, in which the Jacobian is approximated by the following formula:
J (x k ) ⋅ (x k − x k −1 ) ≈ F (x k ) − F (x k −1 )
(5)
The Broyden’s method is computationally more efficient, and its efficiency can be improved by using Sherman-Morison matrix inversion formula [5] for
finding J (⋅) . However, the global convergence properties of this method are not better than those of the Newton-Raphson method. Note that in both methods, the dynamical systems linearize the nonlinear function F ( x ) . This linearization works when the nonlinear nature of the function is not dominant or when successful initial guesses are made. However, the dynamical systems are developed so that their convergence properties remain unknown. Apart from converging to the roots, they may converge to local extrema, diverge, or be unstable, causing numerical overflow. The main aim of this work is to suggest a different approach to the process of finding roots. The synergetic iterative algorithm, which is described in the following section, creates a dynamical system with known properties. Perhaps the single most important feature of this algorithm is that the generated dynamical system will be linear according to the solutions of the system of equations. Therefore, the global convergence of the dynamical system (and of the method itself) is its unique property. −1
3. THE SYNERGETIC CONTROL THEORY APPROACH (SCT) In general, SCT [6,7] provides methods for designing optimal controllers for dynamical systems, where the controllers are coordinated with internal expectations of the systems. The resulting dynamical systems with their controllers have areas of attraction that correspond to the control goals. Introducing such areas of attraction, or attractors, to dynamical systems, is one of the main concepts of the SCT. In this work, the attractors are created at the roots of nonlinear systems of equations. This allows the Synergetic Algorithm (which is described later in this section) to converge and find the roots of the systems. An attractor is a region in the state space of a dynamical system that pulls the trajectories from nearby areas of the state space. Depending on the dimensionality of the state space, attractors can be points, contours, tori or regions of fractal dimensionality. The attractors represent the internal ‘wishes’ of the dynamical system. Whatever the initial conditions are, the system moves towards one of the attractors and remains there infinitely.
If the requirements that the controllers provide to the system cannot be fulfilled, the attractors cannot be created, and the dynamical system becomes unstable. This situation can be identified by stability analysis. If no attractors are present, the system is unstable and it will not converge. This situation can be diagnosed by analyzing stability using the Lyapunov Stability Theory [8, 9] or a potential function of the system [10]. The SCT suggests creating attractors in dynamical systems at those areas of state spaces that correspond to the control purposes. The control purposes are formulated as aggregated macrovariables ψ i that need to be zeroed. The aggregated macrovariables are functions of system variables x and control signals u : ψ i ≡ ψ i ( x , u ) . Hence the attractors need to be introduced where all the aggregated macrovariables are equal to zero. The SCT gives an equation which can be used for creating dynamical systems with attractors at ψ i = 0 .
T ⋅ψ& + ϕ (ψ ) = 0
(6)
where T determines the rate of convergence to the attractor, ψ& is the derivative of the aggregated
macrovariable by time; and ϕ (⋅) - some function that influences the approaching to the attractor. One of possible expressions for ϕ (⋅) is: ϕ (ψ ) = ψ . In this case, the equation is rewritten as:
T ⋅ψ& + ψ = 0
(7)
The solution of this differential equation gives the following function for ψ :
X
0.1
0.05
20
40
60
80
100
Time
-0.05
Figure 2. Convergence of conditions to attractor at ψ
ψ
from different initial
=0
The parameter T in (8), as it can be seen from the equation, determines the rate of convergence. By choosing smaller T , the rate of the transition processes can be increased. Now, let us apply the concepts of the SCT for creating the Synergetic Algorithm for finding roots of the systems of equations. The Synergetic Algorithm Suppose that a system of 1 ≤ j ≤ M nonlinear equations is given as a set of:
f j (x1 , x2 ,K, xN ) = 0
(9)
This set of equations represents the purpose of the control, which is expressed by the aggregated macrovariables. The aggregated macrovariables in this case are equal to the equations: ψ j ≡ f j (⋅) . This assures
(8)
that the attractors will be created at the solutions of the system of equations. The dynamical system is created using equation (7):
where t stands for time. As it is depicted in Figure 2, ψ (t ) is attracted to ψ = 0 from any initial conditions
∂f j (x1,K, xN ) ∂f j (x1,K, xN ) T ⋅ ⋅ x&N + f j (x1,K, xN ) = 0 ⋅ x&1 + (10) ∂xN ∂x1
ψ (t ) = ψ 0 ⋅ e
t − T
ψ0. In order to model the dynamical system, its equations need to be solved for time derivatives. According to time derivatives, these differential equations are linear, thus the solution can be found using Gauss
( ) , or other
Elimination, which has a complexity of O n methods. The resulting dynamical system:
x&i = g i (x1 , x2 ,K, xN , T )
3
(11)
can be iterated using any iterative method, like RungeKutta, etc.
A specific result can be obtained for a one dimensional case:
x& = −
f (x ) T ⋅ f x′ ( x )
(12)
An important property of the dynamical system given by (11), is that its attractors are located where the system of equations has solutions. Besides, although g i (⋅) are nonlinear, the dynamical system is linear according to
ψj.
This guarantees convergence to the
solutions only, provided that they exist. If no roots exist, iterating over the dynamical system is not needed; instead, stability analysis can be used for checking the system’s consistency. An example of such a system is given in Section 4. For cases when there are several roots, several attractors are created in the dynamical system. The initial guess determines the root to which the system converges. By using SCT concepts for solving systems of nonlinear equations, the synergetic algorithm achieves two goals: (1) It checks if the system has roots prior to iterating and (2) It converges to the roots globally, i.e. independently to initial guesses. 4. EXPERIMENTS AND DISCUSSION Examples illustrating the synergetic algorithm are given in this section. The algorithm is discussed in relation to: (1) A simple system of equations. (2) A system that may be inconsistent, (3) A comparison of the the convergence rate versus the Newton-Raphson method. A simple system of equations Suppose that we need to find the roots of the system:
e x + e y = a x− y =b
(
)
1 0.5
Time 20
(14)
Then, the system is solved for x& and y& to obtain the dynamical system equations. Once that is done, the dynamical system is iterated over to obtain the solution.
40
60
80
100
-0.5 -1
Figure 3. Convergence of the system (a= 2, b=1) from different initial values
Figure 3 demonstrates the convergence of the dynamical system for X and Y to the root from several initial guesses. The convergence is finished after t > 60 . The obtained solution matches the root of the system, found analytically: for a = 2 and b = 1 , it is x = 0.38 and y = −0.62 . Inconsistency in equations Let us now discuss the case when the equations have no solution. Without loss of generality and for illustrative purposes only, a one dimensional case is considered here. Suppose that the equation is:
x2 + a = 0
(15)
Depending on the value of a , the equation can have one or two roots, or no roots at all. After applying the procedure to the equation, the dynamical system that converges to the root(s) is:
x& =
(13)
The application of the synergetic algorithm in this case will consist of the following steps. First, the dynamical system based on the equation (10) is constructed:
T ⋅ e x ⋅ x& + e y ⋅ y& + e x + e y − a = 0 T ⋅ ( x& − y& ) + x − y − b = 0
X,Y 1.5
− x2 − a 2 ⋅ x ⋅T
(16)
By introducing a potential function V ( x ) as
∂V ( x ) = − x& ∂x
(17)
we can represent the ‘energy’ of the dynamical system and use it for stability analysis. The set of maxima of the potential function correspond to unstable states (analogous to maxima of ‘potential energy’), minima – to stable states. The potential function for the obtained dynamical system is:
V (x ) =
x2 a + ⋅ ln ( x ) 4 ⋅T 2 ⋅T
(18)
Now consider the three possible situations, based on parameter a . Case 1. a > 0 .
4
V(x)
3
V(x) 2
4 1
2
x
x -4
-2
2
4
-4
-2
2
4
V ( x ) for a = 0
-2
Figure 6. Potential function
-4
When a = 0 , the dynamical system is stable, it has a single attractor and globally converges to the solution as depicted in Figure 7:
-6
Figure 4. Potential function
0.1
V ( x ) when a > 0
0.075
Thus, the potential function has no minima and the system is unstable (Figure 5): X
X
0.05 0.025
Time 2
4
6
8
Time
10
2
4
6
8
10
-0.025
-200
-0.05 -400 -600
Figure 7. Stable system with one attractor when
Case 3. a < 0 . In this case the potential function has two minima as depicted in
-800 -1000 -1200
Figure 5. Unstable system when
a=0
6
a>0
As it is depicted in Figure 5, the choice of close initial guesses causes significant deviations of the trajectories. There is no conversion and the system stops due to numerical overflow. Case 2. a = 0 . In this case, the logarithm is eliminated from the potential function, and the potential function becomes a parabola with a focus at x = 0 as depicted in Figure 6.
V(x)
5 4 3 2 1
x -4
-2
Figure 8. Potential function
2
4
V ( x ) when a < 0
As in the previous case, the dynamical system is globally stable (except for x = 0 ), however now it has two attractors corresponding to the roots. The choice of any of the attractors depends on the initial guesses. In this case, the y axes is a potential boundary between the two attracting regions. A dynamical system cannot intersect this line by itself. Therefore, the choice of the root is determined by the initial guesses, as it is illustrated in Figure 9.
In terms of the SCT, this dynamical system has a bifurcation point at x = 0 , where any deviation from zero causes significant changes in the system. X
4 3 2 1
Time 2
4
6
8
10
-1 -2
This paper discusses the application of the Synergetic Control Theory for solving systems of nonlinear equations. A Synergetic algorithm is presented and its characteristics are discussed. The algorithm has the the following advantages vis-a-vis existing iterative algorithms: (1) The nature of the algorihm is such ahta it makes possible to directly assess the consistency of the systems of equations, (2) The global convergence to the roots of systems is guaranteed and (3) The algorithm exhibits controllable rate of the convergence, which can be made faster/slower than the rate of the other algorithms. BIBLIOGRAPHY
-3
Figure 9. Stable system with two attractors when
a