Polynomial Method for PLL Controller Optimization†

The Phase-Locked Loop (PLL) is a key component of modern electronic communication and control systems. PLL is designed to extract signals from transmission channels. It plays an important role in systems where it is required to estimate the phase of a received signal, such as carrier tracking from global positioning system satellites. In order to robustly provide centimeter-level accuracy, it is crucial for the PLL to estimate the instantaneous phase of an incoming signal which is usually buried in random noise or some type of interference. This paper presents an approach that utilizes the recent development in the semi-definite programming and sum-of-squares field. A Lyapunov function will be searched as the certificate of the pull-in range of the PLL system. Moreover, a polynomial design procedure is proposed to further refine the controller parameters for system response away from the equilibrium point. Several simulation results as well as an experiment result are provided to show the effectiveness of this approach.


Introduction
For nonlinear control systems, one would often like to know the domain of attraction (DoA) of an equilibrium point. Often, this domain is difficult to both find and represent computationally. The usual mathematical tool used for analyzing of the region of attraction is Lyapunov's method. This gives us a sufficient condition for local stability, although it is often difficult to find a Lyapunov function that can be used as a certificate for the whole DoA. Several prior approaches have used quadratic functions, for example [1][2][3]. In particular, the approach of [3] makes use of semidefinite programming to find a quadratic function whose sublevel-set is a good inner approximation to the region of attraction. For system in which the DoA is complicated, an ellipsoid may not provide a good approximation, and the above methods leave a large unexplored region within the DoA.
With recent developments in algebra and sum-of-squares techniques, it is now possible to solve for a Lyapunov function with a more general polynomial form [4,5]. Positive definiteness properties are replaced by sum-of-squares(SOS) constraints which can be efficiently solved using convex optimization. The SOSTOOLS [6] toolbox for MATLAB has been developed as an easy computational tool to solve problems that utilizes the SOS techniques. This approach has also allowed finding a Lyapunov function within some specified semi-algebraic region [7,8]. While Lyapunov approach provides a method to certify a given inner approximation to the domain of attraction, it does not immediately provide a way to find it. Tan [9] later extended this concept by using unions of SOS polynomials to estimate the domain of attraction but Tan's approach needs to solve a bilinear optimization. The level-set method [10] has been developed to find a semi-algebraic representation of the DoA by semidefinite programming. With these polynomial techniques, it is possible to precisely estimate the DoA of a nonlinear polynomial system and to find a suitable Lyapunov function as the stability certificate.
A phase-locked loop (PLL) system is a nonlinear system with limited domain of attraction. Due to its importance in communication systems, analyzing and designing a PLL system has attracted many attentions in this field [11][12][13][14][15][16][17]. The current approach for designing a controller for a PLL system is still largely based on the linear model [17]. Hence, the performance of the resultant system cannot be guaranteed at system states far away from the designed equilibrium point.
In this paper, we utilize the current SOS techniques to analyze the domain-of-attraction of a PLL system. A local Lyapunov function can then be found as the certificate of the DoA using the proposed approach. The Lyapunov function will be further used to improve the stability region and performance of the PLL system. Using this approach, it is possible to design a PLL with a predefined form of controller that has larger domain-of-attraction than the linear design approach. Moreover, since we are designing the system in the nonlinear region, system dynamics outside the linear region can be further refined. Examples of second order PLL systems are used later in this paper to show the effectiveness of this design approach. In the provided examples, we demonstrated that the system designed by the proposed method has 20% larger domain-of-attraction and less overshoot with faster convergent rate than the linear design approach. This paper is organized as follows. Section 2 contains some preliminary knowledge about SOS techniques. The advection algorithms as well as the SOS methods for finding a local Lyapunov function are stated in Section 3. A brief discussion of the configuration of a PLL can be found in Section 4 and the proposed method for analyzing and designing a PLL controller are presented in Section 5. Then the paper is summarized in Section 6.

Preliminaries
The following are some definitions that will be used frequently in this paper. R[x] is used to represent the ring of polynomials in x with real coefficients.
It is also well-known that the converse is not true. Σ denotes the set of all SOS polynomials in R[x]. R + is used to represent the set of nonnegative real numbers. B r is used to represent the open ball with radius r centered at the origin.
Suppose g : R n → R is C 1 . Define the 0-sub-level set of g to be sub(g) ⊂ R n given by Further define the boundary of sub(g) as ∂ sub(g).
One feature of the proposed advection algorithm is that the advection problem can be converted into a semidefinite program. The following is a standard form of a semidefinite program.
where X ∈ R n×n is symmetric. X ≽ 0 means that z T Xz is positive semidefinite for all z ∈ R n .
The condition of one semi-algebraic set containing another semi-algebraic set is one of the key constraints used in this paper. The following lemma shows that this kind of relationship can be converted to constraints on the coefficients of the polynomials. The proof can be found in [4] or [7].
Then sub(q) ⊂ sub(p). Further, given q and the degree bound of p, s 0 , and s 1 , the set of coefficients of p, s 0 and s 1 satisfying (1) is the feasible set of a semidefinite program.
The representation shown in Lemma 1 is one of the simplest cases of Schmüdgen's Theorem [18]. Schmüdgen's Theorem states that if p ∈ R[x] is strictly positive inside a compact semi-algebraic set S generated by p 1 , . . . , p m as S = {x ∈ R n | p i ≥ 0, i = 1, 2, . . . , m}, then , 1} m and s v ∈ Σ. Putinar [19] later showed that under some additional constraints on p i , p has a simpler representation as The gap between Schmüdgen's and Putinar's representation is later investigated by Jacobi and Prestel [20]. In the simple case shown in Lemma 1, if sub(q) ⊂ sub(p) and sub(q) is compact, the representation of p by (1) is always possible. The following result is similar. Given q ∈ R[x], if there exists s 0 , s 1 ∈ Σ and ϵ > 0 such that Usually q is a given polynomial and p is the solution to find such that sub(p) and sub(q) approximately represent the same set with some other constraints on p, such as having lower degree or passing through several pre-specified points. The above results are used to construct such constraints.

Acquiring the Local Lyapunov Function
Finding a local Lyapunov function is coupled with finding the DoA. Without a clear knowledge of the actual shape of the DoA, it is hard to find a Lyapunov function that can be used to represent the entire DoA [4,8]. To deal with this difficulty, we utilize the current development in set advection [10].

Set Advection
In this paper, we will consider the following autonomous systeṁ where f : R n → R n is locally Lipschitz. From the basic local existence and uniqueness theorem [21], given an open subset U ∈ R n , there exist c ∈ R + such that the autonomous system (2) has a unique solution for any z ∈ U in the compact time interval [−c, c]. We define the flow map ϕ t : R n × R → R n to be the local unique solution of For any t ∈ R such that ϕ t (x) exists, the map ϕ t : R n → R n is continuous, invertible and has a continuous inverse, i.e., it is a topological homeomorphism on R n [22].
Given t ∈ R, we define the time t advection operator A t : where C(X, Y ) is the set of functions mapping from X to Y . The map A t is also called the Liouville operator associated with f . A very important property is that it is linear. Figure 1 shows the concept of the advection operator. Given polynomial p, A t maps the coefficients of p to another polynomial q such that sub(q) = ϕ t sub(p). We relate the advection operator to the advection of sets in the following remark.

Time-Stepping
Since we are performing advection, we must use an approximation to the flow map ϕ h with time step h. Many such approximations are possible, and are provided by the theory of numerical integration. The first-order Taylor approximation to where the derivative Dp(x) is a linear map Dp(x) : R n → R n at each point x.
Based on the required accuracy of the advection, we could also choose to use higher order Taylor approximation. However, depending on the system dynamics, this usually will lead to the requirement of using higher degree polynomials in the sum-of-squares constraints. The relationship between the accuracy and the degree of polynomials will be further investigated in future work.

Domain-of-Attraction Estimation
The set advection concept is used to estimate the DoA of a system. We use the following definition of the DoA in this paper.
Definition 1. Suppose f : R n → R n is analytic with the flow map, ϕ, and the origin is asymptotically stable. Define the domain-of-attraction (also called the basin/region of attraction) of f to be R ⊂ R n such that for any x ∈ R, ϕ t (x) is defined for all t ≥ 0 and lim t→∞ ϕ t (x) = 0.
The following properties can be easily derived. The detailed proofs can be found in [10].
Lemma 2. Suppose f is analytic and the origin is asymptotically stable and R ̸ = ∅. Suppose also S 1 ⊂ R and 0 ∈ S 1 , and S 1 is a connected closed positively invariant set. Let h > 0 be a positive constant, and define the backward advection of S 1 to be S 2 , given by Then S 1 ⊂ S 2 ⊂ R, and S 2 is also connected, closed and positively invariant. Furthermore, Theorem 1. Suppose f is analytic and the origin is asymptotically stable and h > 0. Also suppose 0 ∈ S 0 and S 0 ⊂ R is a closed connected positively invariant set, such that there exists ϵ > 0 such that Define the sequence of sets S 0 , S 1 , S 2 , . . . by Then this sequence converges to R in the following sense: (iii) For every x ∈ R, there exists n such that x ∈ S n

Star-Shaped Constraint
For the case of estimating the DoA, we introduce the concept of star-shaped sets. The star-shaped sets have many important properties and can be easily implemented as a semidefinite program. We now start with the first property. The detailed information about the star-shaped set can be found in [10].
Note that a star shaped set S is connected. We now give a simple sufficient condition that ensures a sub-level set is star-shaped. We make the following definition.
Definition 3. Suppose g : R n → R. We call g strictly star-shaped if g is C 1 and further satisfies g(0) < 0 and The following lemma shows the connection between strictly star-shaped functions and star-shaped sets.
Proof. Suppose x ∈ sub(g). Let y : R + → R n be the function The trajectory of y(t) follows the straight line connecting x and the origin. We would like to show that y(t) ∈ sub(g) for all t ≥ 0. We have and since y(0) = x we have g (y(t)) < 0 for all t ≥ 0 as desired. Now we need to show that if y ∈ ∂ sub(g) then there does not exist λ ∈ [0, 1) such that λy ∈ ∂ sub(g). Suppose for the sake of a contradiction that there does exist such a y and λ. We know that there exists such λ > 0, since g(0) < 0. Define the function h : Then the derivative of h is For the purposes of this paper, we would like to construct a convex set of functions whose sub-level sets are connected. Although the convex set of all convex functions on R n will suffice, using it would unnecessarily restrict the class of sets describable to be convex. One cannot simply use the set of all functions whose 0-sub-level set is connected, since this set of functions is not convex. We therefore choose the set of strictly star-shaped functions, which is a convex set. We will use strictly star-shaped polynomials to represent sets. This is significantly more general than existing approaches using quadratic functions [1][2][3]. Also, it has been shown in Lemma 3 that if g is strictly star-shaped, then sub(g) is strictly star-shaped. By using this property, we can easily pose the star-shaped constraints on g to make sub(g) a connected set.

An Algorithm for Backward Advection
Here we will state the result of the backward advection algorithm. Given a strictly star-shaped polynomial g i−1 such that sub(g i−1 ) ⊂ R, and sub(g i−1 ) is bounded and positively invariant, we compute a polynomial g i such that sub(A h g i ) ≈ sub(g i−1 ) as follows.
Pick α > 0 and positive integer d. Solve, using semidefinite programming, the following feasibility problem for Here we introduced an important parameter, α, which we think of as follows. The above algorithm finds a degree d polynomial g i such that g i is strictly star shaped, ϕ h sub(g i ) ⊂ sub(g i−1 ), and ϕ h−α sub(g i ) ⊃ sub(g i−1 ). Hence the parameter α may be thought of as a tolerance that allows for the constraint that g i is required to have degree d or less. Then from the result of Theorem 1, lim i→∞ sub(g i ) converges to the domain-of-attraction. It should be noted that this technique only works in the case that the advected set is positively/negatively invariant.

Stopping Conditions
By using the proposed level-set method, one can successfully propagate the system states backward in time. However, a stopping criterion is still needed to terminate the iterations. To detect the convergence of the advected sets, the closeness of two semi-algebraic sets is analyzed. The following result shows that the closeness of two semi-algebraic sets can be estimated by using scaled sets. The detailed proof can be found in [10].
Theorem 2. Suppose g 1 and g 2 are strictly star-shaped functions, sub(g 1 ) ⊂ sub(g 2 ), and sub(g 2 ) is bounded. Suppose x 1 , x 2 ∈ R n are two points such that and x 2 ∈ ∂ sub(g 2 ), and x 1 = αx 2 for some α ≥ 0. Define the function q : R n → R by Figure 2. Example of stopping conditions. Curve 1 is the inner set and curve 2 is the outer set. Curve 3 is the shrunk set of the outer set. To determine when the algorithm should terminate, one formulates an optimization problem using Lemma 1 to determine the smallest λ > 1 such that sub(q) ⊂ sub(g 2 ), where q(x) = g 2 (λx). Again, this may be evaluated using semidefinite programming. In practice, one picks a λ > 1 in advance, and checks this condition after each iteration. Figure 2 shows an example. Here Curve 1 is ∂ sub(g 1 ), Curve 2 is ∂ sub(g 2 ), and Curve 3 is ∂ sub(q). The largest radial deviation between Curves 1 and 2 is less than 0.3.

The Local Lyapunov Function
We find a local Lyapunov function in order to construct an initial star-shaped positively invariant set. The following result is standard.
Proposition 1. Suppose f : R n → R n is analytic and the origin is a stable equilibrium point. Also suppose V : R n → R is a C 1 function, γ > 0, and the set Then sub(g 0 ) is positively invariant, and sub(g 0 ) ⊂ R.
One simple approach to finding an initial sub-level set is to find a quadratic Lyapunov function for the linear model of the system, and use a small sub-level set of this quadratic polynomial as the initial set.
An alternative method which often gives a much larger initial set is as follows. Choose a polynomial p ∈ R[x] such that sub(p) ⊂ R. We then solve the following convex feasibility problem. Find V ∈ R[x] and s 0 , s 1 ∈ Σ such that Similar methods for finding local Lyapunov functions along with details on the construction of the associated semidefinite program may be found in [5,8]. Here we have added the first constraint to ensure that V − γ is strictly star-shaped for γ > 0. Note that these constraints imply that all sub-level sets of V are compact. Given V , we then solve the convex program where ϵ > 0 is small. The optimal γ satisfies sub(V −γ) ⊂ sub(p). Then V and γ satisfy the assumptions of Proposition 1 and so we may use g 0 = V − γ as the function defining our initial level-set. After we have found an estimate of the DoA, we can then use Proposition 1 to find a Lyapunov function that can be used to describe the behavior of the system within the DoA. To adequately describe the system behavior, this approach requires that we have a good estimate of the DoA. The following examples show that the level-set algorithm can precisely estimate the DoA.

Examples of DoA Estimation
Example 1. Consider the following dynamical systeṁ The origin is a locally stable equilibrium point. Here we start with the initial polynomial g 0 = 2x 2 + 2y 2 − 1. The results of the level-set method are shown in Figure 3, using time step h = 0.2. It can be seen that the successive iterates approach the true DoA.
This system is locally stable around the origin. An initial sub-level set given by the quadratic polynomial g 0 = 4x 2 + 4y 2 − 1 is used which can be verified to be positively invariant. A time step of h = 0.2 is used. The time tolerance parameter α is 0.02. The even-numbered iterates g 0 , g 2 , g 4 , . . . are shown in Figure 4. For reference, some of the iterates are listed below and normalized to allow integer coefficients. It can be seen that the iterates gradually approach the exact boundary of the DoA. After 30 iterations, the solution covers most of the stable region.
After 40 iterations, the stopping criteria allowing an absolute radial change of 0.01 has been met. The final result is shown in Figure 4 as Curve 4. Curve 1 is ∂C(g 0 ) if we were using the semidefinite-programming based procedure in Section 3.7. For comparison, Curve 2 is the result of [2] and Curve 3 is the result of [1]. Example 3. The following dynamic system is the Example S4 taken from [3].
Here the initial set is a small circle around the origin. After 20 iterations, the estimated boundary of the DoA has reached the pre-specified bound. The final iterate is shown in Figure 5. The dashed curves show the two system trajectories which are used to represent the true boundary of the domain of attraction. Curve 2 represents the result of the final iterate. Curve 1 is the result from [3].   Figure 6 shows the basic configuration of a PLL. It has three components; a phase detector, a loop filter, and a voltage controlled oscillator(VCO). The VCO generates an output signal whose phase, θ 0 (t), depends on the phase, θ i (t), of the input signal. The PLL is phase locked when the phase error ϕ(t) = θ i (t) − θ 0 (t) is a constant value and the loop is in stable equilibrium state. Usually, it is desired that the phase error, ϕ(t), is maintained at zero. Of interest is the behavior of the phase error ϕ(t). Because of its sinusoidal nonlinearity in the PLL, the phenomenon of chaos is believed to exist [11,12] and its inherent chaotic behavior for broadening the pull-in range of PLL has also been realized [13,14]. A nonlinear controller can drive PLL from chaotic state into periodic state or vice versa [15]. For higher-order PLL, it is not possible to determine whether the loop will or will not slip cycles using the initial frequency alone. In this case, one might define the pull-in range as the separatrix ordinate at ϕ = 0 [17]. Analyzing the DoA of the PLL system provides a better description of the region in which a PLL locks up without slipping. The Lyapunov method has been used for stability analysis in control systems. Here the advection algorithm will be used to find the guaranteed stability boundary of the PLL system and the associated local Lyapunov function is then used to further refine the controller parameters. In [16], a Lyapunov styled analysis for PLL system up to third order is presented. The method shown in this section provides a way to analyze the DoA for a more general system. Also, the form of the Lyapunov function used here is much more flexible. Figure 7 shows the nonlinear model of the PLL. The sine function here represents the phase detector of the system. K in Figure 7 stands for the loop gain of the system. F (s) is equivalent to the low pass filter shown in Figure 6 and it corresponds to the controller of the PLL. Finally, the integrator in Figure 7 is the voltage or numerically controlled oscillator. The key idea of a PLL system is to use the command, y 2 , from F (s) to steer the oscillator such that θ 0 (t) tracks θ i (t) as closely and quickly as possible.

Second Order PLL
To use the nonlinear design approach, start with a reference design. The reference design used in this paper is the linear model of a PLL system. A Proportional-Integrator (PI) controller is chosen to be the filter, F (s), as Using the model shown in Figure 7, it is routine to check that the resulting dynamic equation of the system is Assume that the received signal frequency is varying linearly with time and has zero radial acceleration and let x 1 = ϕ, x 2 =φ. The PLL system can be rewritten as the following state space model Equation (5) is the final nonlinear model of the second order PLL. A linearized model can then be derived asẋ The filter F can then be designed using existing linear design approach [17]. One typical choice is to let ω n = 15, ζ = 0.707, K = 1, where ω n is the natural frequency, ζ is the damping ratio, and K is the overall gain. The two coefficients of the filter are then given as τ 1 = K ω 2 n , τ 2 = 2ζ ωn

Phase-Locked Loop Analysis and Design
In this section, the second order PLL controller will be used to demonstrate the nonlinear design approach. The same approach can also be applied to the third order controller design.

Pull-In Range of the Traditional PLL System
Now the advection algorithm can be applied to the PLL nonlinear system. To reduce the required number of iterations, a local Lyapunov function is used as the initial set. After a few iterations of the algorithm, it gives us the estimated domain-of-attraction of the system. The result is shown in Figure 8. Note that the estimated region is based on the Taylor series expansion of the sine and cosine functions. From Ston-Weierstrass Theorem, we can approximate the sine and cosine functions to the desired accuracy within a bounded interval. In this example, two degree-10 polynomials are used to approximate sine and cosine functions and the estimated region is only valid between −π and π. More terms of the Taylor series could also be used to improve the accuracy.

PLL System Controller Design
After getting a good estimated DoA, a local Lyapunov function that describes the system behavior can be easily computed using Proposition 1. This local Lyapunov function can also be used to describe the trajectories of the system in the DoA. Figure 9 shows several Lyapunov level-sets obtained from the SOS approach.  The SOS techniques can be applied to the design of a controller. For the PLL system, it is desired to design a system which has a larger DoA or faster converging speed. Here, the objective is to find possible system parameters such that the same Lyapunov function is still valid and the system converges faster. This can be done by solving a semidefinite program.
Suppose V is a local Lyapunov function for the PLL system. Using the SOS technique, solve the following optimization problem: where a, b ∈ R + specify the domain of the constraints and q is a positive definite performance polynomial specified by the user. p(k) is a linear constraint of controller parameter k. As before, s 1 , s 2 , s 3 , s 4 are SOS polynomials.
Since V is now a given function, the above constraints are linear in the controller parameters. The first constraint shows that V is a valid Lyapunov function in sub(V − a). This constraint is used to specify the desired DoA to maintain. The second constraint along with the objective function will put an upper bound on the derivative of the Lyapunov function in sub(V − b). Faster decreasing rate implies faster converging speed. This specifies the performance requirement of our system. The user could also put different performance requirements in different sub-level sets of V .
Besides dynamic performance constraints, noise bandwidth constraints will be applied as well. The noise bandwidth for this PLL system with PI controller has the following form [17] Assume ζ ≥ 1, ω n ≥ 1. Then This linear upper bound will be used to find a set of controller parameters that have better dynamic performance while maintaining the same noise bandwidth. The system phase portrait as well as the estimated stable region are shown in Figure 10. The nonlinear design has K = 1, ω n = 10.813, and ζ = 1.3303. The noise bandwidth is 8.1082 Hz, which is slightly higher than the noise bandwidth of the linear design, 7.9546 Hz. From the phase portrait, the nonlinear design has approximately 20% larger guaranteed domain-of-attraction. It can also be observed from the phase portrait that the nonlinear design has less overshoot than the linear design. This shows that the proposed method increases the system performance while not sacrificing too much of the noise rejection capability.
A Simulink model is used to compare the nonlinear designed PLL controller with the linear design. In this Simulink model, the sinusoidal input is collapsed by measurement noise and clock noise with zero mean and variances 0.1 and 0.0001, respectively. Figure 11 is the phase error of the two designs. It is clear that the nonlinear designed controller has a much faster convergence rate than the original design. Both systems are tested on a NORDNAV R25 software GPS receiver. This receiver collects, down-converts and samples the GPS data by the front end, so that the collected GPS data can be post-processed repeatedly using different tracking-loop filter orders. The results are shown in Figure 12.
It can be seen that the original system has some overshoot and converges around 300 ms. The nonlinear design has much less overshoot and converges about five times faster than the original design.

Summary
In this paper, we presented a method of designing a PI controller of a PLL system. This design approach is based on the polynomial nonlinear model of the PLL system. This approach starts with the linear design of the controller and then estimates the DoA of the linear designed system to get the suitable local Lyapunov function for the system. The Lyapunov function is then used as the performance constraints to further refine the performance of the system outside the linear region. The DoA of the initial design can also be extended to get a better pull-in region of the PLL system. This approach gives us a way to design a fixed form controller for a nonlinear system.