VIEWS: 12 PAGES: 43 POSTED ON: 8/19/2011
Analog and Digital Control 1 Lecture 3 Kirsten Mølgaard Nielsen Based on lecture notes by Jesper Sandberg Thomsen 1 Outline Discrete equivalents Overview (emulation – finding a discrete equivalent to the continuous system) Numerical integration Zero-Pole matching (mapping) Hold equivalents Design by emulation Evaluation of design Nyquist stability criterion Continuous case Discrete case 2 Discrete Equivalents - Overview r(t) e(t) controller u(t) plant y(t) + D(s) G(s) - Translation to discrete plant Translation to discrete controller (emulation) Zero order hold (ZOH) Numerical Integration • Forward rectangular rule • Backward rectangular rule • Trapeziod rule (Tustin’s method, bilinear transformation) • Bilinear with prewarping Zero-Pole Matching Emulation Hold Equivalents Purpose: Find a discrete transfer function which • Zero order hold (ZOH) approximately has the same characteristics over • Triangle hold (FOH) the frequency range of interest. Digital implementation: Control part constant between samples. Plant is not constant between samples. 3 Numerical Integration Fundamental concept Represent H(s) as a differential equation. Derive an approximate difference equation. We will use the following example Notice, by partial expansion of a transfer function this example covers all real poles. Example Transfer function Differential equation U ( s) a H ( s) u au ae E ( s) sa 4 Numerical Integration TransferFunction U (s) a H (s) E (s) sa Differential Equation u au ae u (t ) au( ) ae( ) d t (continous) 0 kT T au aed kT T au aed kT u (kT ) 0 u (kT T ) area of au ae over kT T , kT 5 Numerical Integration Now, three simple ways to approximate the area. kT-T kT Forward rectangle approx. by looking forward from kT-T Backward rectangle approx. by looking backward kT-T kT from kT Trapezoid approx. by average kT-T kT 6 Numerical Integration Forward rectangular rule (Euler’s rule) (Approximation kT-T) Difference Equation u1 (kT ) u1 (kT T ) T au1 ( kT T ) ae( kT T ) (1 aT )u1 (kT T ) aTe(kT T ) Transfer Function U1 ( z ) (1 aT ) z 1U1 ( z ) aTz1E ( z ) U1 ( z ) aTz1 aT a H F ( z) E ( z ) 1 (1 aT ) z 1 z 1 aT z 1 a T 7 Numerical Integration Backward rectangular rule (app kT) Difference Equation u2 ( kT ) u2 ( kT T ) T au2 ( kT ) ae( kT ) u2 ( kT ) aTu2 (kT ) u2 ( kT T ) aTe( kT ) 1 aT u2 ( kT ) u2 ( kT T ) e( kT ) 1 aT 1 aT Transfer F unction 1 aT U2 ( z) U2 ( z) E ( z) (1 aT ) z 1 aT U2 ( z) 1 aT aTz a H B ( z) E ( z) 1 1 (1 aT ) (1 aT ) z 1 z 1 a (1 aT ) z Tz 8 Numerical Integration Trapezoid rule (Tustin’s Method, bilinear trans.) (app ½(old value + new value)) Difference Equation u3 ( kT ) u3 ( kT T ) au3 ( kT T ) ae( kT T ) au3 ( kT ) ae( kT ) T 2 1 ( aT / 2) u3 ( kT T ) aT / 2 e(kT T ) e(kT ) 1 ( aT / 2) 1 ( aT / 2) Transfer F unction U3 ( z) aT ( z 1) a HT ( z) E ( z ) ( 2 aT ) z aT 2 2 ( z 1) a T ( z 1) 9 Numerical Integration Comparison with H(s) H ( s) Method Transfer function a a Forward rule H F ( z) sa z 1 a T a a Backward rule H B ( z) sa z 1 a Tz a a Trapezoid rule HT ( z) sa 2 ( z 1) a T ( z 1) 10 Numerical Integration Transform s ↔ z Method Approximation z 1 Forwardrule s or z 1 Ts T z 1 1 Backward rule s or z Tz 1 Ts 2 ( z 1) 1 (T / 2) s Trapezoid rule s or z T ( z 1) 1 (T / 2) s Comparison with respect to stability In the s-plane, s = jw is the boundary between stability and instability. 11 Numerical Integration - stability Forward z 1 Ts for s jw 1 Backward z for s jw 1 Ts 1 (T / 2) s Trapezoid z for s jw 1 (T / 2) s (Forward) (Trapezoid) (Backward) Region of Region of Region of s-LHP s-LHP s-LHP 12 Numerical Integration Some conclusions Higher distortion with the forward and backward rule. The forward rectangular rule could cause a stable continuous filter to be mapped into an unstable digital filter. The trapezoid rule maps the stable region in the s- plane exactly into the stable region of the z-plane. Though, the jw-axis is compressed into the unit circle with length 2p. The trapezoid rule appears superior. But, let us make some further investigations..... let us calculate the power |H(jw|2 at, for example, w = a – the change in frequency response. 13 Numerical Integration a H ( jw ) jw a a2 1 H ( jw ) 2 2 2 2 w a 2 (w / a ) 1 (w a ) 1 H ( ja) 2 2 What about HT(z) ? at which frequency is |HT(z)|2 = 0.5 ? 14 Numerical Integration a 1 HT(z) as HT ( z) 2 1 2 before but re- j tan(wT / 2) a j tan(wT / 2) 1 T aT formulated with z e jwT . Power wh en z1 e jw1T (page 194) 2 12 1 H T ( z1 ) 2 1 2 2 tan(w1T / 2) 12 a T 2 2 aT a tan(w1T / 2) or w1 arctan( ) T T 2 I.e. high 2 aT sampling Notice, w1 ( ) a , aT / 2 1 frequency T 2 15 Numerical Integration Trapezoid with pre-warping maintain the power at some specified frequency 1. Choose a critical frequency w1 2. Replace w1 by a´ such that new a = a’, 3. and insert in H(s). a’ = (2/T) tan(w1 T/2) 4. Find the discrete filter H(z) based on the modified continuous transfer function H(s). s H ( z) H ( ) w1 s 1 z 1 w1 tan(w1T / 2) z 1 w w1 H ( z) H ( j1) 16 Numerical Integration Plots for T = 0.1 ws ≈ 60 wb Example, 3th order filter H (s) 1 s 3 2s 2 2s 1 with 0.1 T 1 2 Prewarpingat w / wb 1 17 Numerical Integration Plots for T = 1 ws ≈ 6 wb Notice, prewarp 18 Numerical Integration Plots for T = 2 ws ≈ 3 wb Notice, prewarp 19 Zero-Pole Matching Exact map between poles (z = esT) exists. Basic idea, simply use same map for zeros. Rules 1. All poles are mapped by z = esT. For example, s = -a maps to z = e-aT 2. All finite zeros are mapped by z = esT. For example, s = -b maps to z = e-bT 3. Basically, zeros at jw = infinity maps to z = ejp = -1 (representing the highest frequency) 4. Identical gain at some critical freq. (typically, s=0) H(s) at s=0 = H(z) at z=1 20 Zero-Pole Matching Example (1) Compute the discrete equivalent by zero-pole matching a H(s) sa Pole (rule 1) s a z e aT Zero (rule 3) H ( s ) s 0 z 1 z 1 a 1 e aT Gain (rule 4) K 1 K z e aT z 1 s a s 0 2 z 1 H ( z) K z 21 Zero-Pole Matching a zp-match z 1 H(s) H ( z) K sa z s-plane z-plane 22 Zero-Pole Matching Including delay in computer Notice, u(k) depends on e(k). This comes from rule 3 U ( z) z 1 1 z 1 H ( z) K K E( z) z 1 z 1 u(k ) u(k 1) Ke (k ) Ke (k 1) Solution, first zero for s=infinity is mapped to z=infinity (as in Matlab). This gives U ( z) 1 z 1 K K E( z) z 1 z 1 u(k ) u(k 1) Ke (k 1) 23 Hold Equivalents Purpose, design a discrete filter H(z), such that input e(k) is samples of e(t) H(z) has an output u(k) u(k) approximates the output of the continuous filter H(s) whose input is the continuous signal e(t). Zero order hold (ZOH) (step invariant), approach below e(t) e(k) eh(t) u(t) u(k) Sampler Hold H(s) Sampler 24 Hold Equivalents Hold operation q(t ) 1(t ) 1(t T ) 1 s Q ( s ) 1 e Ts E h ( s ) Q ( s ) e( k ) e( k ) is a gain Q(s) e(t) e(k) eh(t) u(t) u(k) Sampler Hold H(s) Sampler 25 Hold Equivalents ZOH transformation Notice, we can find H(z) by e(t ) (t ) e(k ) k Take inverse Laplace to find a continuous signal. Sample this continuous signal to find a discrete signal. U ( s) H ( s) Q( s) e(k ) H ( s) Q( s) (k ) u(t ) L H ( s) Q( s ) (k ) L H ( s) 1 e Ts 1 1 1 s (k ) e(t) e(k) eh(t) u(t) u(k) Sampler Hold H(s) Sampler 26 Hold Equivalents Continuous output signal 1 1 s u(t ) L H ( s) 1 e Ts (k ) 1(t) 1(t T ) L1 Hs(s) (k ) Discrete output signal and transfer function H(z) u(k ) u(t )kT t H ( s) Z u(k ) U ( z ) H ( z ) (1 z 1 ) Z s H ( s) 1 H ( s ) where Z Z L s s 27 Hold Equivalents Triangle hold equivalent Extrapolate the samples to connect sample to sample in a straight line. 1.0 Impulse response for triangle hold e(k ) u (t ) Hold System T T eTs 2 e Ts H (s ) Ts 2 Result ( z 1) H ( s ) H ( z) Z 2 Tz s 28 Hold Equivalents Same third order filter as before + Zero-pole o ZOH x Triangle T=1 (next T=2) 29 Hold Equivalents Same filter as before + Zero-pole o ZOH x Triangle T=2 30 Hold Equivalents Concluding remarks Methods, we have looked at Numerical integration methods Zero-pole matching (mapping) Hold equivalents All methods works fine for high sample frequency Matlab 31 Evaluation of Design Discrete controller by emulation Can be calculated using one of several techniques. Evaluation discrete equivalent (design) We need to analyze in the z-domain. Process expressed in z-domain (ZOH) Notice, other eval. might be necessary (simulation) 32 Evaluation of Design Example, using Matlab Cont. sys.(antenna) Cont. controller 1 (10s 1) G( s) D( s) s(10 s 1) s 1 sysD_c = tf([10 1],[1 1]); T = 0.2; % around 30*wn. Also use, T = 1 sysD_d = c2d(sysD_c,T,’matched’); % zero-pole matched sysG_c = tf([1],[10 1 0]); sysG_d = c2d(sysG_c,T,’zoh’); % zero order hold sysCL_c = feedback(sysD_c*sysG_c,1); sysCL_d = feedback(sysD_d*sysG_d,1); step(sysCL_c); hold on; step(sysCL_d); 33 Evaluation of Design T = 0.2 (ws ≈ 30 wn) 34 Evaluation of Design T=1 (ws ≈ 6 wn) The sample delay gives phase lag. Thus, less PM which gives higher overshoot 35 Nyquist Stability Criterion Check closed loop stability with controller KD(s) and plant G(s) Y ( s) K D( s ) G ( s ) b( s) , D( s) G( s) R ( s ) 1 K D( s ) G ( s ) a( s) Characteristic equation (CE) a( s) Kb( s) 1 K D( s) G( s) 0 0 a( s) Notice, number of poles in CE equals number of open- loop poles in D(s)G(S). (known) Notice, zeros in CE are closed-loop poles! (unknown) 36 Nyquist Stability Criterion Approach, check characteristic equation Let us follow a contour encircling the RHP. Basic idea: Evaluate CE = 1+KD(s)G(s) at the contour, and draw image. (We can just as well evaluate KD(s)G(s) and look at -1 ) K D( s) G( s) Notice, poles or zeros gives encirclements 37 Nyquist Stability Criterion Criterion (continuous case) Z=P+N Z = number of unstable closed-loop poles. Thus, for stability we want Z=0. (unknown) P = number of unstable open-loop poles of KDG (known). N = total number of encirclements of the point -1 by evaluation of KDG (taken in same direction as s. Thus, N might be negative). 38 Nyquist Stability Criterion Interpretations and explanations If the s-plane contour encircles a zero of 1+KDG in a certain direction, the image contour will encircle the origin in the same direction. (related to Z) If the s-plane contour encircles a pole of 1+KDG in a certain direction, the image contour will encircle the origin in the opposite direction. (related to P) The net number of same-direction encirclements, N, equals the difference N = Z- P. Actually, we should investigate 1+KDG and encirclements around the origin, but easier to investigate KDG and encirclements around -1. 39 Nyquist Stability Criterion Discrete case Same approach, almost... Unstable region is outside the unit circle. Easier to evaluate using the stable region. 40 Nyquist Stability Criterion Discrete case CE: 0 = 1+KD(z)G(z) Total number of poles = n Number of stable zeros of 1+KD(z)G(z) is n-Z Number of stable poles of 1+KD(z)G(z) is n-P Criterion (discrete case) Z = (n-Z) – (n-P) or Z=P–N 41 Nyquist Stability Criterion Discrete case Determin the number P of unstable poles of KDG Plot KD(z)G(z) for the unit circle. This is the counter- clockwise path around the unit circle. Set N equal to the number of Counter-clockwise encirclements of the point -1 on the plot Compute Z=P-N. The system is stable if and only if Z=0 42 Nyquist Stability Criterion Example (7.9) Evaluate the stability of the unity feedback discrete system with the plant G(s) with sampling time T=2 and ZOH. K=1 Use controller KD(z) = K. 1 G( s ) s( s 1) Solution (matlab) sysC = tf([1],[1 1 0]); sysD = c2d(sysC,2); nyquist(sysD); The system is openloop stable and there are no encirkelments of -1 => the system is closed loop stable for K=1 43