2.3 Modified First Order Simulation Methods 

2.3.1 Simplified Solutions Based on Approximations of e^{z} 

If t_{N} and t_{D} are constant (and assuming a constant Dt), then the recurrence coefficients A and B are also constant and can be precomputed offline using the optimum constantt formulas given by equations (2.2.34). However, if t_{N} and t_{D} are variable, the optimum coefficients given by (2.2.44) and (2.2.45) are rather cumbersome, and require storage of past values of t_{D} and t_{N}. We have shown that the variablet response can be approximated by the constantt solution if Dt is small enough, but in this case the A and B coefficients must be recomputed online for each interval, and it may be desirable to simplify the computation. Since the most timeconsuming aspect of the computation is the exponential in equation (2.2.34), the most common approach is to substitute some approximation for this exponential. 

If Dt/t_{D} is much less than 1, we can substitute the first terms of the Taylor series expansion for e^{z} into equation (2.2.34), which leads to the family of approximate recurrence coefficients presented in Table 1. 

Table 1 

Many lead/lag algorithms correspond to one of these approximate solutions, or some variation. For example, one simulation program uses an algorithm equivalent to 

_{} 

This corresponds to the lowest order approximation from Table 1, except that B includes a Dt in the numerator. The effect of this Dt is to introduce a bias in the response to nonconstant input functions (because B is the coefficient of Dx in the recurrence formula). This bias was not evident during testing because the only evaluation of the algorithm was based on its response to step inputs, during which Dx was always zero except at a step. 

Another example is a commercial embedded software package that uses a lead/lag algorithm with the coefficients 

_{} 

Here, the solution corresponds to the second order approximation from Table 1, except for the –Dt/2t_{D} component in the right hand term of B. 

It’s interesting to compare all five of the modified B coefficients discussed so far for a particular case. Consider a lead/lag with t_{N}/t_{D} = 10 and Dt/t_{D} = r < 1. The three levels of approximation from Table 1 are 

_{} 

whereas the dynamic simulation algorithm of the first example gives B = 10 + r, and the embedded commercial algorithm of the second example gives B = 10 – (9/2)r – (1/4)r^{2}. According to the Taylor series approximation method, the former approximation is wrong in the first order correction for nonvanishing r, while the latter approximation is wrong in the second order. The reason that the designers of these algorithms chose to deviate from the Taylor series approximation as shown here is not known. 

Instead of the Taylor series approximation, we could make use of the series expansion of the natural log (see Section 2.2.5). If we take just the first term of this series and make the substitution 

_{} 

we have 
_{} 

Taking the exponential of both sides gives 

_{} 

In the range of interest, this is a much more efficient approximation of the exponential than the truncated Taylor series. Substituting this for the exponential in equations 
(2.2.34) gives the recurrence coefficients 

_{} 

Incidentally, some commercial applications have implemented simple firstorder lags which compute y_{c} based only on x_{c} and y_{p}. Since they do not store the past value of x, they effectively assume x_{p} = x_{c}. With this assumption the value of B is unimportant, because it is the coefficient of x_{c} – x_{p} which vanishes. Hence we’re only interested in the A coefficient, which these applications take to be A = Dt/(t_{D} + Dt). Clearly the optimum representation of a firstorder lag to this level of approximation would be with A = Dt/(t_{D} + Dt/2), as can be seen from the preceding formula. It is not known why the application designers chose to neglect the factor of 1/2 on Dt in the denominator. This makes the effective time constant of the algorithm equal to t_{D} + Dt/2. Typically Dt is about 1/4 of t_{D}, so this represents up to 25% error in the time constant. 


2.3.2 Finite Differences and Bilinear Substitution 

The most direct approach to determining an approximate recurrence formula for a differential equation of the form 

_{} 

is to substitute simple finite differences for the differentials about the nearest convenient point. The derivative of the function y(t) at the point halfway between the past and current instants can be approximated by 

_{} 

The value of the function y(t) at this point is approximately 

_{} 

Similarly for the function x(t) at the same instant we have the approximations 

_{} 

Making these substitutions into equation (2.3.21) and rearranging gives the recurrence formula 

_{} 

By equations (2.2.23), the standard form recurrence coefficients are 

_{} 

These are identical to equations (2.3.12), which shows that the direct application of finite differences gives the same firstorder recurrence formula as results from substituting the exponential approximation of equation (2.3.11) into equations (2.2.34). 

There is another approach that, called bilinear substitution (or Tustin’s algorithm), that leads to the same result as finite differences (and the optimum exponential substitution). To apply this method we first rewrite equation (2.3.21) in the form 

_{} 

where D is the differential operator. We then make the substitution 

_{} 

where the symbol z^{1} is the “past value” operator. For example, z^{1} y_{c} = y_{p}. Substituting equation (2.3.25) into equation (2.3.24), we have 

_{} 

Multiplying through by Dt(1+z^{1}) gives 

_{} 

Rearranging terms, we have 

_{} 

Carrying out the indicated past value operations gives 

_{} 

Solving for y_{c}, we arrive at equation (2.3.22). Thus, when applied to a firstorder transfer function, the bilinear substitution method gives the same recurrence formula as the method of finite differences. However, this is not the case for the second order functions, as is shown in Section 3.3. 


2.3.3 PreWarping 

A technique called “prewarping” is sometimes applied to lead/lag algorithms in an attempt to correct the amplitude of the output responding to an oscillating input signal. The standard recurrence coefficients for this method are 

_{} 

Recalling that as z goes to zero, tan(z) approaches z, we see that the prewarped coefficients reduce to equations (2.3.23) if Dt is small relative to t_{N} and t_{D}. 

The rationale for prewarping is clarified by the following example. Given a lead/lag simulation with timestep size Dt = 0.25 seconds and a sinusoidal input function with a period of T = 1 second, each cycle of the input will be described by four discrete values, as illustrated in Figure 1. 



Clearly a sample will rarely be taken at an absolute peak of the wave, even if T is not an exact multiple of Dt and the samples precess through the cycle. Therefore, assuming linear interpolation between the discrete values, we will truncate most of the peaks and underrepresent the true amplitude of the input function, which results in a corresponding underestimation of the output amplitude. In an attempt to correct for this, prewarping biases the recurrence coefficients to reduce the effective lag of the filter, thereby increasing its output amplitude for a given oscillating input. However, the bias necessary to cancel out the underrepresentation of the input varies with the ration of T to Dt, so the prewarped coefficients are really only correct for one particular frequency – by convention, the “half power frequency”. At lower frequencies, and for most nonoscillatory inputs, the prewarping coefficients yield less accurate results than the approximate coefficients given by equations (2.3.23). Hence, prewarping is a highly specialized technique, intended only for a limited range of frequencyresponse applications. 


2.3.4 Composite Past Value Recurrence 

All the solution methods discussed so far have assumed the availability of xp and yp, that is, past values of the input and out functions. This requires the dedication of two RAM memory locations for each lead/lag in the system. In contrast, the composite past value (CPV) recurrence scheme required only one stored value. In the constantt case this recurrence scheme is functionally equivalent to the standard recurrence formula for any given solution coefficients A and B. However, when CPV recurrence is applied in the variablet case, it does not give exactly the same results as the standard formula. 

To derive CPV recurrence formula, we begin by rewriting the standard formula as follows 

_{} 

Dividing and multiplying the two pastvalue terms by 1B gives 

_{} 

If we define the quantity in the brackets as the “composite past value”, denoted by z_{p}, then we can rewrite equation (2.3.41) as 

_{} 
or 
_{} 

The composite current value is defined as 

_{} 

To simplify this, we substitute the expression for the y_{c} from equation (2.3.42), which gives 

_{} 

Combining terms in x_{c} and simplifying, we have 

_{} 

Therefore, the CPV recurrence formulas are 

_{} 

The derivation of these formulas does not distinguish between the past and current values of A and B, thereby tacitly assuming that these coefficients are constant. However, when t_{N} and t_{D} are variable, the standard coefficients A and B are also variable, in which case the CPV recurrence formulas are not exactly equivalent to the standard formula on which they are based. To demonstrate this, we can solve equation (82.3.43) for z_{p} and substitute this into equation (2.3.44) to give 

_{} 

A and B are now subscripted with a “c” to denote the current values, since only the current values are available each time the recurrence formulas are executed. Equation (2.3.45) is also valid at the past instant, which is to say 

_{} 

If we substitute this expression for z_{p} into equation (2.3.43), the resulting equation for y_{c} can be reduced to the form 

_{} 

where A_{eff} and B_{eff} are the effective standard recurrence coefficients, related to the actual coefficients used in the CPV recurrence formula by 

_{} 

Because of this effect, there are circumstances in which CPV recurrence yields unrealistic output response, even if the optimum formulas for A and B are used. The simplest illustration of this is the case where t_{N} and t_{D} initially have different values and then become equal. The true first order response is for the output to asymptotically converge on the input function once t_{N} = t_{D}, and this can be accurately simulated using the standard recurrence formula. However, this response cannot be reproduced using CPV recurrence. To understand why, notice that whenever t_{N} = t_{D} the current value of B (denoted by B_{c}) must equal 1, as is clearly shown by equation (2.2.34). Substituting this into equations (2.3.46) gives the effective coefficients A_{eff} = 1 and B_{eff} = 1. Inserting these into the standard recurrence formula gives the result 

_{} 

Thus, using CPV recurrence, the output undergoes a “step” change and “instantaneously” becomes equal to the current input, regardless of the output’s past value. 

A possible objection to equation (2.3.41) is that it appears to predict divergent behavior when B_{p} equals 1, whereas in practice CPV recurrence formulas never yield divergent response. The explanation is that B_{p} can only equal 1 if B_{c} of the previous iteration was 1. But we have already seen that if B_{c} was 1 for the previous iteration, then y_{c} was unconditionally set equal to x_{c} during that iteration. Therefore, for the present iteration, y_{p} must equal x_{p}. Recalling that A is the coefficient (x_{p} – y_{p}) in the standard recurrence formula, we see that, in effect, the divergence of A_{eff} at B_{p} = 1 is always “cancelled” by the vanishing of (x_{p} – y_{p}). 


2.3.5 Discussion 

Digital simulations of continuous dynamic functions have a constraint on their execution frequency imposed by the intrinsic information rate of the input and output signals. The execution frequency must be high enough to resolve and convey the significant characteristics of these signals. However, the simplified methods of simulation described in Sections 2.3.1 and 2.3.2 have another limitation: The approximations on which they are based are valid only if the execution period Dt is small relative to the characteristic times of the function (i.e., t_{N} and t_{D}). To illustrate that these two limitations are logically independent, consider the simple lag response 

_{} 

where t = 0.05 seconds. Equations (2.3.23) give the recurrence formula 

_{} 

Now suppose that, based on their maximum information rates, we have determined that x and y can be adequately represented with a sampling frequency of 4 Hz, corresponding to a Dt of 0.25 seconds. Substituting this into equation (2.3.51) gives 

_{} 

A sample of the response given by this recurrence formula for a single timestep is illustrated in Figure 2. 



Figure 2 

The output computed by equation (2.3.52) crosses over the input, which is obviously not representative of a pure lag response. In general, recurrence formulas based on equations (2.3.23) will not give an accurate representation of a lead/lag function unless Dt is less than 2t. This means that the execution frequency of the present example would have to be raised from 4 Hz to at least 10 Hz. In contrast, if the optimum coefficients given by equations (2.2.34) are used, the recurrence formula is 

_{} 

for which an execution frequency of 4 Hz is adequate. Thus the use of nonoptimal coefficients imposes unnecessarily high execution frequency on the simulation. 

Another way of looking at this is to regard the error in the computed value of y_{c} for any given interval Dt as consisting of two parts: one part due to the uncertainty in the input function, and one part resulting from the approximations made in determining the coefficients of the recurrence formula. The first part is unavoidable (for a given Dt), but he second part can be eliminated entirely by using the optimum recurrence coefficients. 

With the technique of prewarping, an attempt is made to use these two errors to compensate for each other. However, as discussed in Section 2.3.3, these two errors offset each other only under a limited range of conditions. Essentially, prewarping attempts to account for nonlinearity in the input signal over the interval Dt, but there is really no general way of determining the nonlinearity of a signal based on just two samples (x_{p} and x_{c}). In applying its “correction”, prewarping makes use of information about the input signal beyond what can be inferred from the two given points. Specifically it depends on the knowledge that the input is periodic with a given frequency. Without such special information, a general method of accounting for nonlinearity of the input signal must make use of at least three values of x. With three values we can determine the optimum recurrence formula based on second order interpolation of the input function. Denoting the “past past value” of x by x_{pp}, the expanded recurrence formula can be written 

_{} 

where 
_{} 

By the same reasoning that led to the relation f_{1} + f_{2} + f_{3} = 1, we have the relation 

_{} 

Using this, we can rewrite equation (2.3.53) in a generalized standard form analogous to equation (2.2.22) as follows 

_{} 

where 
_{} 

If x_{pp}, x_{p}, and x_{c} are colinear vs time, then the quantity multiplied by C vanishes, and equation (2.3.54) reduces to equation (2.2.22). 

Regarding the technique of CPV recursion, it is worth noting that the form developed in section 3.3 is not unique. Beginning with equation (2.3.40), we could have defined the composite past value directly as just the sum of the past value terms 

_{} 

Then the CPV recurrence formulas would have been 

_{} 

The advantages of equations (2.3.43) and (2.3.44) is that they can be executed with just two multiplications, whereas equations (2.3.55) and (2.3.56) involve three multiplications. This is significant because multiplication is a timeconsuming operation in most digital devices. However, equations (2.3.55) and (2.3.56) will give better accuracy in the variablet case because they are equivalent to using the effective standard coefficients 

_{} 

These coefficients never yield singular values of the coefficients, unlike equations 
(2.3.43) and (2.3.44) when B_{c}=1, so no special provisions are needed to account for this computational singularity (recognizing that it cancels when multiplied by (x_{c}  x_{p})). 
