Dynamic Zero Finding for Algebraic Equations

In a variety of contexts, for example the solution of differential games and the control of power systems, the design of feedback control laws requires the solution of nonlinear algebraic equations: obtaining such solutions is often not trivial. Motivated by such situations we consider systems of nonlinear algebraic equations and propose a method for obtaining their solutions. In particular, a dynamical system is introduced and (locally) stabilizing control laws which ensure that elements of the state converge to a solution of the algebraic equations are given. Illustrative numerical examples are provided. In addition it is shown that the proposed method is applicable to determine the equilibria of electrical networks with constant power loads.

H 2 /H ∞ control problems have been considered.General algebraic equations have been considered in [17].
In this paper we consider general nonlinear algebraic equations for which we seek a solution.To this end a dynamical system, which is such that at any of its equilibria components of the state coincide with a solution of the algebraic equations, is introduced.Feedback control laws which, subject to certain conditions, locally asymptotically stabilize an equilibrium of the system are provided.The solution of the underlying algebraic equation can then be determined from the steady-state values of the trajectories of the closed-loop system.The main results of this paper can be summarised as follows.The problem of solving a system of algebraic equations is posed as a feedback stabilization problem for a dynamical system.A systematic method of constructing (locally) asymptotically stabilizing control laws is then proposed.The latter contribution draws inspiration from [18].Therein the authors provide a method for constructing stabilizing control laws for systems with a particular structure with applications to power systems (see also [6], [19]).More precisely, the authors consider dynamical systems with a triangular structure and provide stabilizing feedback control laws through the introduction of an augmented system.The results presented in this paper are applicable to a wider class of problems.
The remainder of this paper is organised as follows.The problem considered is introduced in Section II.The main result, that is the design of a control law rendering the equilibria of a closed-loop system (locally) convergent to a solution of this problem, is provided in Section III.Remarks and numerical examples concerning some particular classes of algebraic equations are then provided in Sections IV and V.The results are then illustrated on a power systems problem in Section VI before some concluding remarks are given in Section VII.
Notation: The sets of natural numbers, real numbers and complex numbers are denoted by Z, R and C, respectively.The empty set is denoted by ∅.I n denotes the n×n identity matrix.Given x ∈ R n we define |x|2 = x T x, and for a matrix A, A denotes its transpose, vec(A) denotes its vectorization and adj(A) denotes its adjoint matrix.Given a square matrix A, σ(A) denotes its spectrum.The Kronecker product between two matrices A and B is denoted by A ⊗ B. All functions are assumed sufficiently smooth.The gradient of a function f is denoted by ∂f ∂x .The Euclidean basis vectors are denoted by e i ∈ R n , i = 1, . . ., n.The relations vec(ABC) = (C ⊗ A)vec(B) , where A, B, and C denote matrices of appropriate dimensions, are used repeatedly.

II. PROBLEM FORMULATION
Let f : R n → R n denote a function described by1 where z ∈ R n , q ∈ R n and P : R n → R n×n .Assume there exists z * ∈ R n such that and Suppose, furthermore, that the image of f around z * is onto a neighborhood of 0. Note that (3) implies that z * is an isolated solution of (2).Define the dynamical system where x(t) ∈ R n , z(t) ∈ R n and Θ(t) ∈ R n×n are the states of the system and u(t) ∈ R n and W (t) ∈ R n×n denote control inputs.For this system we consider the problem of finding a static state feedback such that the closed-loop system has an asymptotically stable equilibrium at (0, z * , P (z * )).Thus, the second component of the steady-state of any trajectories converging to this equilibrium provide a solution to the algebraic equations f (z) = 0.

Remark 1:
The assumption that the image of f is onto a neighbourhood of 0 around z * expresses nothing else than Brockett's necessary condition for the feedback stabilization, via smooth feedback, of system (4), see [20].

A. Design of the control law
To provide a solution to the aforementioned problem we introduce a variable z d ∈ R n , which is implicitly defined by the relation where φ : R n → R n is selected such that For instance one could select φ(x) = tanh(x), for > 0, or φ(x) = x.Note that, provided Θ(t) is invertible for all t ≥ 0, The main contribution of this paper, provided in Section III-C, is the construction of a feedback control law which locally asymptotically stabilizes an equilibrium of the system (4)- (5).A Lyapunov-like design approach is adopted to construct the control law.To this end, define the state vector X := col(x, z, vec(Θ)) ∈ R 2n+n 2 and consider the function where The time derivatives of each of the component functions of V along the trajectories of the system (4)-( 6) are given by V1 where we have defined the n × n 2 matrix and the n 2 × n matrix and used the fact that żd Now, define the matrix It is clear from (10) that selecting the control u and vec(W ) to satisfy where k z > 0 and k Θ > 0 are tuning parameters, yields Note that if we select φ(x) such that x φ(x) It is clear that the derivations above rely on the assumptions that Θ and M (X) are invertible matricesan issue that is addressed in the following two subsections.

Remark 2:
The assumption that Θ is invertible for all t ≥ 0 can be interpreted as an observability condition on the output z d implicitly defined in (5).

B. A preliminary lemma
In this subsection we give a verifiable condition for the invertibility of M (X).Towards this end, define the desired closed-loop equilibrium where z * ∈ R n is a solution of (2), the sets and the mapping . .e n adj(Θ)(q + δ)adj(Θ)]P z (z) , (16) where δ ∈ R n is a parameter.
The following assumption is made to ensure the (local) existence of u and W satisfying (11).
Assumption 1: There exists > 0 and µ > 0 such that Assumption 1 and the following lemma enables the design of the locally asymptotically stabilizing control laws provided in Section III-C (and in Section IV and Section V for linear and quadratic problems, respectively).

C. A stabilizing controller
To begin with we make the important observation that, in view of condition (3), the equilibrium X e is in the set N .This allows to design a controller which locally stabilizes the equilibrium X e as established in the following proposition in which denotes the level sets of the function V (X) defined in (8), with z d given in (7), and k = max{k : Thus, Vk denotes the largest level set of V contained in the set Ω µ .
Proposition 1: Consider the dynamical system (4)-( 5) and a function φ(x) satisfying (6).Suppose the condition (3) and Assumption 1 hold and let X(0) ∈ Vk.Then the dynamical system (4)-( 5) in closed-loop with the control law u and W satisfying (11) is such that the equilibrium X e given in ( 13) is locally asymptotically stable and, in particular, the following implication is true Remark 3: In [18] a method for designing stabilizing controllers for a class of triangular systems has been proposed.The method therein considers systems similar to (4) in which P (z) is diagonal.The result provided in Proposition 1 is applicable to such systems, as well as a wider class of nonlinear systems.
Two special cases, namely systems of linear equations and systems of quadratic equations are considered in the following sections.Whereas in the former case P (z) is a constant matrix, in the latter P (z) is linear in z.Numerical examples are provided for both cases, which provide useful insights into the results expressed by Proposition 1.

IV. SYSTEMS OF LINEAR EQUATIONS
Consider the case in which f defined in (1) is linear.In this case, P (z) = P is a constant matrix that, in view of (3), is full rank.We seek to obtain a solution of the system of linear equations Recall the notation introduced previously, namely X = col(x, z, vec(Θ)) and X e = col(0, z * , vec(P (z * ))).Corollary 1: Consider the system (4)-( 5) with P (z) = P and a function φ(x) satisfying (6).Let Θ(0) be such that det(Θ(0)) det( P ) > 0 .
Then there exists u and W satisfying (11) and the system (4)-( 5) in closed-loop with u and W is such that X e is an asymptotically stable equilibrium and lim t→∞ z(t) = z * .
Remark 4: Corollary 1 establishes that for the case in which P (z) is a constant matrix, the Assumption 1 in Proposition 1 is trivially satisfied provided det(Θ(0)) and det(P (z * )) have the same sign, see (18).

Remark 5:
In Proposition 1 it is assumed that P is full rank, such that (17) has a unique solution.If instead P ∈ R m×n , with m < n, ( 17) is an underdetermined system and may have no solutions or infinitely many solutions.In the latter case Corollary 1 may be generalised to achieve local asymptotic stability for an equilibrium corresponding to a solution of the underdetermined system.

A. A numerical example
A system of linear equations, as studied in Section IV, is considered and Corollary 1 is used to obtain the solution of (17).
Considering the case in which n = 20, the parameters P and q are randomly generated in MATLAB 3 such that P is full rank.The system (4)-( 5) in closed-loop with the feedback control law satisfying (11) is simulated with 4 φ(x) = x, k z = 1 and k Θ = 1 and with the initial condition z(0) = 0, z d (0) = 0, Θ(0) = I n and x(0) = q.Note that the initial condition is such that det(Θ(0)) det( P ) > 0 as required in ( 18) and Corollary 1.The time histories of the components of x are shown in Figure 1 and the time histories of the components of z (top) and z d (bottom) are shown in Figure 2.
V. SYSTEMS OF QUADRATIC EQUATIONS Consider now the case in which (1) is quadratic, i.e. the case in which P (z) defined in ( 9) is linear in z.In this case P z (z) = Pz is a constant matrix.

Θ(0)
z d (0) Corollary 2: Consider the case in which ( 1) is a system of quadratic equations.Let z * denote a solution of (1) such that the condition (3) is satisfied.Suppose for all |δ| < , and for some > 0. It follows that Assumption 1 is satisfied.
Remark 6: For the case in which ( 1) is quadratic, since P z is constant, Assumption 1 is a condition on Θ(0) (and P (z * )).

A. A numerical example
Consider the quadratic equation where q 1 ∈ R and q 2 ∈ R. Consider the selection P (z) = −z 2 0 z 1 z 2 , and consider the case in which q 1 = 8 and q 2 = 16.In this case the system of equations (20) has two solutions.Following the main result given in Proposition 1, the system (4)-( 5) in closed-loop with the feedback given by ( 11) is simulated for different selections of parameters.Let Θ 0,1 = diag{−1, 1} and Θ 0,2 = diag{1, −1} and let z d0,i = [0, i], for i ∈ Z.The initial condition Θ(0) is varied between Θ 0,1 and Θ 0,2 , whereas z d (0) = z d0,i , for i = −4, . . ., 4. The considered initial conditions are referred to according to Table I.The remainder of the parameters are selected as follows 5 The trajectories of z are shown in Figures 3 for the initial conditions I 4 (solid, black line), I 3 (solid, dark grey line), I 2 (solid, light gray line), I −4 (dashed, black line), I −3 (dashed, dark grey line) and I −2 (dashed, light grey line).The solid, square markers denote the initial conditions, whereas the solutions of (2) to which the trajectories converge are denoted by the diamondshaped markers.

VI. APPLICATION TO ELECTRICAL NETWORKS WITH CONSTANT POWER LOADS
In this section we apply the developed numerical method to compute the steady state of an m-port linear time-invariant AC electrical network with Constant Power Loads (CPLs) working in sinusoidal steady-state at a frequency ω 0 ∈ R + , see [5] and references therein for additional details on this problem.
The description of this regime in the frequency domain is given by the relation where 6 V ∈ C m and I ∈ C m are the vectors of generalized Fourier transforms of the port voltages and currents, respectively, and K ∈ C m captures the effect of the external (current or voltage) AC sources, all evaluated at the frequency ω 0 .Defining the apparent power at the CPLs as S i := P i + jQ i , i ∈ {1, . . .m}, (22) 6 To simplify the notation the argument jω 0 is omitted in the sequel.
where P i ∈ R and Q i ∈ R are the active and reactive power at port i, respectively, the CPLs constraint the network via Vi where (•) c denotes the complex conjugate operator.Obviously, a necessary and sufficient condition for the existence of a sinusoidal steady-state (at a given frequency ω 0 ) is that the complex equations (21), ( 22) and ( 23) have a solution.It has been shown in [5] that this set of equations can be written in the form q − P (z)z = 0 , where q := col(P, Q), z := col(Re{I}, Im{I}) and with K = col(Re{K}, Im{K}) and Hence, the problem of finding the equilibrium current I can be solve using the method presented in this paper.

TABLE I THE
SELECTION OF INITIAL CONDITIONS CONSIDERED IN THE SIMULATIONS.