The effects of proportional steering strategies on the behavior of controlled clocks

Two approaches to steering clocks to a reference standard via frequency steers have been explored in the literature, the pole placement (PP) approach and the linear quadratic Gaussian (LQG) approach. The PP method sets the response time of the controlled clock. LQG minimizes a user-selected cost function, which is a linear combination of the expected variances of the phase, frequency, and control effort (frequency steers). This paper unifies these techniques by developing an analysis that quantifies what time constants follow from the LQG-derived gains, and what costs are incurred by the PP-derived gains. With this analysis, one can compute how sensitive the outcomes would be to variations of the steering constants. Another outcome we present is a quantification of the degradation due to a common form of suboptimal state estimation, as well as that due to delays before steering decisions can be implemented into the system.


Introduction
Control theory is often taught at the graduate level, appears in generalized form in textbooks and publications such as [1][2][3][4], and has been applied to disciplining oscillators [5,6]. The theory is not yet complete because it does not fully address the effects of suboptimal modelling, small variations in the loop gains, and other properties [7,8], although simulations can provide a useful supplement in this regard [9]. As will be shown, the commonly-made mistake of assuming the measurement is noiseless results in steering following measurement noise, and this oversteering will degrade the performance. Another outcome of poorly chosen gains might be unexpected oscillatory behavior.
This paper extends published analyses of the PP and LQG techniques [5], which are capable of finding the gains that optimize performance as defined by either setting the time constants (PP) or by minimizing a linear combination of the phase, frequency, and control variances respectively (LQG). We unify the two techniques by deriving the formulas for quantitative estimates for the variances and time constants when the gains are at the desired value, as well as for all other values that do not result in unbounded performance. Unfortunately, it is not always possible to correctly or even optimally estimate the state, and in some systems there may be long delays before a steer can be implemented. We therefore derive formulas to quantify the effects of suboptimal state estimation when measurement noise is ignored, and for data processing lags.
We follow the Kalman formalism (appendix A) for state estimation since it yields an optimal estimate of the current state under very general assumptions, such as independent noise effects and proper parameterization.
Although a variety of control schemes are possible, we consider only those in which the clock is controlled by shifting its frequency as fractions (gains) of its phase and frequency deviations from the reference standard; typical examples are a laboratory timescale steered to coordinated universal time (UTC), or a precise clock steered to a local time scale via a microstepper. A justification for ignoring the frequency drift is that frequent enough steers will naturally remove it with little impact. Section 1 briefly describes the LQG approach to finding the gains that produce minimum cost. Section 2 shows what the resulting cost variances are for every gain pair that leads to stable behavior. Section 3 is the PP approach and shows how the poles of the Z-transform are related to the time constants and performance variances of the system. Section 4 quantifies how suboptimal state estimation degrades the steering. Section 5 shows how to quantify the effects of a time-lag in the post-measurement steering decisions. Appendix A provides a review of the Kalman formalism, and appendix B shows how to compute the expected Allan Variance of the controlled clock.

Optimization via linear quadratic Gaussian (LQG) approach
The LQG approach finds the value of the gain vector G that minimizes any desired linear combination of the phase, frequency, and frequency-steer variances [1][2][3][4][5]. The coefficients of those variances are termed costs and written W Q (1,1), W Q (2,2), and W R respectively. The phase cost has dimension 1/τ 2 , while the other two costs are dimensionless. Of the three cost parameters only two are independent, since an overall scaling of all three parameters does not affect the minima. A reason to minimize the control variance is because the controlled clock is usually more stable than the reference for small τ , and excessive steering will destabilize short-term performance.
As a consequence of the well-known Separation Theorem, which states that the optimal control strategy is independent of the optimal state estimation method, optimal gains are independent of the measurement and process noise, although they do depend on the sampling interval. As shown in the references, the formalism to compute the gain results in a Riccati equation whose solution can found using the Matlab function dlqr. The LQG solutions can therefore be summarized in figures 1 and 2, which show the phase and frequency gains as a function of the costs, scaled so that the control cost is 1. While there is no dependence on an overall scaling of costs, there is an implicit dependence on the unit interval through the cost definitions. In subsequent sections we compute and present the costs as a function of the gains, and it will be seen that the indicated minima are consistent with the LQG-derived ones in figures 1 and 2.

The expected variances of the phase, frequency, and frequency steers as functions of the gain vector, in steady-state
In this section we first define the terms for the Kalman filter ( [1][2][3], and appendix A). Following [5], we assume a system whose 'true' and estimated states x true,k and x k are stable and defined by the difference between their phase p k and frequency f k with respect to a reference clock. After the measurement at time t k , the estimated states are x k = Ç p k f k å , and their covariances are P k . The steering goal is p k = f k = 0. The measurements are at times t k , spaced by intervals τ. The vector x −,k is the Kalman filter's estimate of the true state function x true,k at an instant just before the measurement at t k , when its covariances are given in the matrix P −k . The measurement is assumed to be immediately used to create the state estimate x k , and the frequency steering, based on that estimate, is immediately implemented. We utilize standard Kalman formulas such as Here H is the 2 × 1 measurement vector, which for our purposes is (1,0). K is the Kalman gain, ν k and w k are the measurement noise at time t k and the process noise leading to time t k respectively, ∆x k = x true,k − (x − ) k , and z k is the innovation. The Kalman gain is computed using the process noise variance Q and the measurement noise variance R, and the state  Using the Kalman formalism, it is possible to compute the expected variances of the phase, frequency, and steers for any set of gains that yield finite steady-state results in the long term. We first find the steady-state Kalman gain K ∞ , and the steady state pre-and post-measurement covariances P −∞ and P ∞ , by requiring them all to be unchanged through a measurement cycle and applying appendix A's equations (A.8), (A.10) and (A.19). They lead to a discrete-time algebraic Riccati equation for P −∞, which can be solved using, for example, the Matlab call [P −∞ , L, F] = dare Φ T , H T , Q, R . After enough iterations to reach steady-state, the state evolution equation can be written with all the gain-dependent terms absorbed into This is a Lyapunov equation, in Matlab: Since the noise parameter ν k and the prediction error ∆x k are uncorrelated with each other and (x − ) k (though not x k ): (2.9) The variance of the control, u, is The derivation of the expected Allan variances (and Hadamard variances) is a straightforward extension of the formulas developed in this section, and provided in appendix B. For gains large enough, or equal to 0, the Lyapunov equation does not yield a solution. In these cases, one or more of the variances   is unbounded; the formulas delineating the region of stability are derived in the next section. By way of illustration, figure 3 shows the RMS as a function of gains for a process noise of 0.1 and measurement noise of 0.1; normalized so its minimum is 1. In the LQG formalism, this plot would correspond to a cost function weighting only the phase variance. The minimum occurs at a phase and frequency gain of 1, and the cost goes up more rapidly as the gains are varied further from the optimum. Figure 4 shows that the frequency RMS would be reduced by 70% if the phase gain is lowered to near 0 (with a phase gain of exactly 0 the phase cost would be unbounded). Figure 5 shows the cost defined by control effort (RMS of steers) will strongly decrease as the gains approach 0. More complicated cost functions result in different behavior, as the total cost is a linear combination of the phase, frequency, and control variances.
As an example, figures 6-8 show the phase, frequency, and steering of a controlled clock whose intrinsic frequency varies with process noise Q = 0.01, and whose measurement noise variance is R = 0.01. For both gain vectors plotted, (0.2, 1.0) and (0.2, 0.3), the mean phase and frequency of the steered clock is zero.

Pole placement (PP)
The essence of PP is to choose the time constants by appropriately placing the poles in the Z-transform. One first decides what time constants are desired, and these correspond to the poles. The desired value of the poles then yields the gains. In order to develop this approach, we explicitly work with the phase and frequency components of the state vector (p k and f k ), and note that f k is (p k − p k−1 )/τ. X(z), the Z-transform of a time series, p k , can be written: For exponential decay with time constant T, Including the effects of steering in the two-state phase and frequency model: where the frequency y k−1 = (p k -1 − p k−2 )/τ and n k is the combined effects of the uncorrelated process noise and optimally-filtered measurement noise from k = 0. Taking the Z-transform of both sides: (3.8) The 'extras' incorporate the initial conditions and the noise, while the poles and their associated time constants T pole are given by We note that equation (3.9) can also be derived from the eigenvalue equation Ax k = e −τ /Tpole x k , and this is developed in section 5. As in [5], the solution is critically damped (T = T c ) if the discriminant of (3.9) is zero, and then If the discriminant is negative there will be two oscillation frequencies. Expressing the pole components as re iθ = e −τ /T e iθ = e −τ /T e 2πifoscτ , where f osc is the frequency of the oscillations, the decay constant T is given by the magnitude of (3.9):

12)
and the frequency is f osc If the discriminant is positive, there are two time constants and the magnitude of their associated terms will depend upon the details of the disturbances that excite them. If one of the poles is negative, there is oscillatory behavior at the Nyquist frequency f osc = 1 2τ . Poles near the origin of the unit circle (|z pole | = 1) have a smaller value of e − τ T , and therefore a faster response and a smaller decay time T. A pole on the origin of the unit circle has an 'infinitely fast' response. Poles near the unit circle have a slower response (large time constant T), consistent with the  gains being low. Poles outside the unit circle are unstable. Lower frequency gains yield oscillatory behavior while larger gains have two associated time constants. Figure 9 shows how the positive, negative, mixed, and imaginary pole configurations are a function of phase and frequency gain. The unit circle is mapped onto the diagonal line segment which marks the boundary of stability, and the line segments along the gain = 0 axes from the origin to the diagonal line segment. The term 'critically damped' is generally applied only when the gains are <1, however we apply it to all cases wherein the discriminant is zero, even if the pole is negative and therefore oscillatory.
The fact that the pole must have an absolute value <1 can be used to explicitly delineate the region of stability in figure 9, defined as where a finite impulse yields a finite response. The largest absolute value of z pole occurs when the two terms being summed on the right-hand side of equation (3.9) have the same sign. If the first term is negative and the discriminant is positive, the equation at the point of instability reduces to (3.15) Squaring both sides, and assuming the discriminant is positive, If the first term in equation (3.9) is positive and the discriminant is positive and the boundary is at g 1 = 0.
If the discriminant is negative, squaring the real and imaginary parts of equation (3.9) shows that the boundary is at g 2 = 0. Figure 10 shows the frequencies, in units of 1/τ, when the poles are imaginary. These are exponentially damped with a time constant that depends only on g 2 , as given in equation (3.12).
Using the formulas of section 3, it is possible to plot the interplay between the critical gains and their associated time constants (figure 11), as well as the RMS associated with them ( figure 12). The figures illustrate the generally greater influence of frequency gains, and that a long time constant   minimizes control effort at the price of very large phase RMS, while a short time constant degrades the frequency performance.

Suboptimal estimates
The Separation Theorem only implies that the optimal gains derived by LQG are optimal when the state itself is optimally estimated. An obvious example would be a state measured with a system that always returns the same value no matter what the state is; in that case setting both gains to zero may be the best no matter what performance is desired. A more common error is to assume that the measurements of the state are noiseless, and that the frequency is exactly found through differencing the last two measurements. It is easily shown that this is equivalent to setting the Kalman gain K = Following the formalism already developed in section 2, one can replace the steady state Kalman gain with the suboptimal gain K and follow all the equations before (but not including) (2.9) to obtain: (4.1) In this case the predicted state is no longer uncorrelated with the prediction error ∆x k , and two terms that were zero in the derivation of equation (2.9) are no longer so: where Σ x−∆x = (x − ) k ∆x T k and Σ ∆x = ∆x k ∆x T k . The reason the second and fourth terms in (4.2) are not zero is because, while the optimal prediction is just as likely to differ positively as negatively from the true value no matter what its sign is, our suboptimal estimate's difference with the true value is more likely to have the same sign. Using equation (2.4), and later (2.7): Equation (4.7) is a Sylvester's equation, solvable in Matlab as Assembling the now known components, one can now solve equation (4.3) as a Lyapunov equation: As before, the control effort is given by u 2 = GΣ x G T .
(4.10) Unfortunately, the above formalism is invalid for the special case where g 2 = 1, because the matrix A = is singular. However, there is only one unknown and the equations can be directly solved. The ratios shown in figures 13-15 are the RMS for the suboptimal case divided by the optimal RMS, for two different cases of process and measurement noise. In our 2-state form ulation, the fractional degradation depends on the ratio of the measurement noise to the process noise. The consequences of the noiseless approximation become more serious as the ratio R/Q increases. For R = 0.1 and Q = 1, all ratios are within 5% of optimal.
In some cases, the ratios are less than unity for low gains, particularly for the phase RMS. This is partly because the low gains in effect average over the noise much as an optimal state estimator would filter them out, and also because the mis-estimation of the frequency and phase would be correlated for every point, which leads to larger steers and so the gains are effectively larger-thus reducing the phase and frequency offsets faster.

Steering when there is an information lag
In steering a master clock to UTC, the information needed to steer comes in 30 d batches, approximately 10 d late. A similar situation exists for Rapid UTC (UTCr), which is computed in 7 d batches, arriving about 3 d late. We approximate this by setting τ equal to the batch size, and assuming a steer is implemented only after a fractional lag λ < 1. Any available measurements between times k and k − 1 can be used to estimate x k , so long as they take into account the steers s k+λ applied within that interval. The prediction of the state at the time the information is received (and implemented) is based upon propagating the last known state forward to the time the information is delivered (with Φ λ ) and the steers can be denoted s k+λ 3) The phase and frequency variances are derived as in equation (2.9), if one replaces A by A λ , and the control variance is In order to study the associated time constants, it is possible to develop the eigenvalue approach, which is complementary to the Z-transform method provided in section 3. Ignoring measurement and process noise, and denoting the eigenvalues as E: The determinant equation is The solution to this is identical to equation (3.9), with the eigenvalues E = e −τ /Tpole . Since this expression does not depend on the lag λ, the time constants would be the same function of the gains as derived previously. Figures 16-18 show the ratio of phase, frequency, and control RMS, as a function of gain, to their values when there is no lag. The ratios are still closer to unity for clocks that have lower process noise. While reducing the rather modest degradations due to the phase lag are one reason for the large improvements observed when laboratories steer to UTC using UTCr, the principle reason is that the weekly spacing of the UTCr deliveries entails considerably less process noise between estimates than is associated with the monthly spacing of the UTC deliveries via the Circular T.
Note that while the phase costs for non-zero lags are larger than for zero lag for gains <1, they are not always larger otherwise. This can be understood as follows: if the phase and frequency offsets have the same (positive) sign, the steer will be more negative due to the phase increase over the lag, thus tending to counteract the frequency more and lowering the contribution of the next measured point to the phase variance. It is equally likely that only the frequency offset is negative and then the computed steer would be less, tending to increase the contribution of the next measured point to the variance. The symmetry is broken because the variance depends on the square of the offsets, and it is therefore more sensitive to changes in the larger values that happen when phase and frequency offsets have the same sign. Even so, the reduction in variance is only with regards to the measurement points; inserting the steers in between measurements adds instabilities within the intervals between measurements.
When data are steered to UTC or UTCr, the data come in batches and this analysis would be strictly correct if only the last datum of each batch were used. In such cases, the user may wish to use the intermediate points to estimate the final state, and apply either the LQG or PP approach to determine what the phase and frequency of the state should be at the time when the next batch is slated to arrive. The actual steers could be those that would bring about that final state with the minimum amount of control, as in [6] and [10]. This has been termed 'gentle steering'.

Conclusions
Although the need for simulations always remains, particularly for non-Gaussian behavior, the analytic techniques needed for a complementary analysis of the pole placement and linear quadratic approaches have been demonstrated. These allow estimates of the quantitative degradation due to suboptimal state estimation and delays in information transfer, such as the latency inherent in the monthly computations of UTC and weekly computations of UTCr.

Disclaimer
The software package Matlab is mentioned for technical clarity and because other workers may find references to it to be useful. As a matter of policy, neither the US Naval Observatory nor the US Government will endorse any commercial product.

Appendix A. The Kalman formalism
This section is a review of the Kalman filter [1][2][3], in order to define the notation in the presence of proportional steering in frequency. Following [5], we assume a system whose 'true' and estimated states x true,k and x k are stable and defined by the difference between their phase and frequency with respect to a reference clock, x k = Ç p k f k å , and their covariances are P k The steering goal is p k = f k = 0. The discrete time series of state estimates, measured at times t k , separated by intervals τ, is denoted x k . The vector x −,k is the Kalman filter's estimate of the true state function x true,k at an instant just before the measurement at t k , when its covariances are given in the matrix P − , k . The measurement is assumed to be immediately incorporated into the state estimate x k , and the frequency steering is based upon that estimate, using the standard Kalman formula Where H is the 2 × 1 measurement vector, K is the Kalman gain, ν k is the measurement noise, and ∆x k = x true,k − (x − ) k . The units are arbitrary in the sense that the 'time coordinate' can be turns of phase, nanoseconds, or radians. The frequency and control (steering) terms have units of the time coordinate divided by the sampling interval τ, which is taken to be 1 in the figures.
For measurement k, Immediately before and just after incorporating the measurement at estimation time k, the steering, state estimate, and its covariance are as follows:  The Allan Variance of the frequency steers is found by multiplying the gain vector times the difference between x k and equation (A.3)'s expression for x k+n with the transpose of that product: All of these expressions can be easily though tediously coded on a computer. As a check, the computed Allan variance should approach the frequency variance divided by τ, for large τ.