Tag: Frequency Response

Discretization Techniques

In Digital Filter-Implementation we have shown how digital filters are actually implemented using difference equations. In the blogs, Filters and Controllers in control system, we have shown a couple of basic filters and the PID controller in the continuous time domain. Finally, in Digital Control Systems 1.01 and 1.02 we shown what comes along with digital systems.

In this article we show how to transfer, transfer functions defined in the continuous-time domain to discrete-time domain; discretization. There are several techniques to accomplish this; i) modeling the sampled system or ii) by using Z-transform mappings.

  1. Modeling the sampled system
    • Zero-order hold
    • First-order hold
    • Other higher order hold circuits
  2. Z-transform mappings
    • Forward euler
    • Backward euler
    • Trapeziod / Tustin
    • Tustin with pre-warp
    • Matched pole zero
    • Modified matched pole zero

In this article we only go in detail in the Z-transform mapping techniques.

Finite difference methods

The first three methods actually use the finite difference methods to approximation the derivative and obtain a mapping from the continuous (s = j\omega_c) to the discrete (z = e^{j\omega_dh}) time domain. Herein, the forward and backward Euler methods use the following approximations

    \[\dot{x}(kh) \approx \frac{x(kh +h) + x(kh)}{h} \rightarrow s = \frac{z - 1}{h}\]

    \[\dot{x}(kh) \approx \frac{x(kh) - x(kh - h)}{h} \rightarrow s = \frac{1 - z^{-1}}{h}\]

Whereas the trapeziod method, also referred as the Tustin’s method, uses the approximation

    \[\frac{\dot{x}(kh) + \dot{x}(kh + h)}{2} \approx \frac{x(kh + 1) - x(kh)}{h} \rightarrow s = \frac{2}{h}\frac{z-1}{z+1}\]

Note that the forward Euler, stable poles in s may be mapped to unstable poles in z. For both backward Euler as Tustin’s method, stable poles in s implicates stable poles in z.

Tustin with pre-warping

When determining discrete-time frequency response of the Tustin rule, H_d(z) with z = e^{j\omega_dh}, to a continuous time transfer function H_c(s), it reveals a relation between the continuous and discrete time domain; Every frequency point in continuous time domain is mapped to the discrete time domain via \omega _{c}={\frac{2}{h}}\text{tan}\left(\omega_dh/2\right) \Leftrightarrow \omega_d = \frac{2}{h}\text{arctan}\left(\omega_ch/2\right).

    \begin{gather*}H_d(z) = H_c(s)|_{s = \frac{h}{2}\frac{z-1}{z+1}} \Leftrightarrow H_d(z) = H_c\left(\frac{h}{2}\frac{z-1}{z+1}\right) \\ \Leftrightarrow H_d(e^{j\omega_dh}) = H_c\left(\frac{h}{2}\frac{e^{j\omega_dh}-1}{e^{j\omega_dh}+1}\right) \Leftrightarrow H_d(e^{j\omega_dh}) = H_c\left(j\frac{2}{h}\text{tan}\left(\frac{\omega_dh}{2}\right)\right)\end{gather*}

The tangent introduces a non-linear mapping where the distortion is small when \omega_d h is small. We can compensate for this by setting \omega_c =\frac{2}{h}\text{tan}\left(\omega_dh / 2\right) for every frequency we have control over. That is to say, first replace the critical frequencies in the continuous-time transfer functions via

    \begin{gather*}(s + \omega) \rightarrow (s + \hat{\omega}) \text{ with } \hat{\omega} = \frac{2}{h}\text{tan}\left(\omega h / 2\right) \\(s^2 + 2\beta\omega s + \omega^2) \rightarrow (s^2 + 2\beta\hat{\omega} s + \hat{\omega}^2) \text{ with } \hat{\omega} = \frac{2}{h}\text{tan}\left(\omega h / 2\right)\end{gather*}

Secondly apply the normal Tustin transformation s = \frac{2}{h}\frac{z-1}{z+1}. Finally, scale the static gain of the discrete transfer function to match the gain of the continuous-time transfer function, e.g., H_c(j\omega_0) = k \cdot H_d\left(e^{j\omega_0 h}\right) with \omega_0 = 0 and solving for k.

Remark that most programs, such as Matlab only offer to replace one critical frequency via their discretization (c2d) function. They apply the transformation

    \[s = \frac{\omega_0}{\text{tan}\left(\omega_0 h / 2\right)}⁡ \frac{z - 1}{z + 1}\]

In which \omega_0 is the frequency chosen by the user to be matched in both continuous as discrete time domain.

Matched Pole Zero

This method maps the poles and zeroes according the the relation z = e^{sh}. That is to say,

    \begin{gather*}(s + \omega) \rightarrow z - e^{-\omega h} \\(s + a)^2 + b^2 \rightarrow z^2 - 2\left(e^{-ah}\text{cos}\left( b h \right)\right)z + e^{-2ah}\end{gather*}

If the numerator (zeros) is of lower order than the denominator (poles), add powers of (z + 1) to the numerator, until numerator and denominator are of equal order. This causes an averaging of the current and past input values, as in Tustin’s method.

However, it might be possible that it is desirable to have the output only depend on the past values. For instance, due to costly computing power. In that case, only add as much powers of (z + 1) such that the numerator is of lower order than the denominator by 1. In which case, this method is referred to as the modified matched pole zero.

Likewise, as with the Tustin transform with pre-warping, scale the static gain of the discrete transfer function. Note that this method preserves stability and provides less high frequency response error.

Frequency response example

The figure below shows the frequency response of a continuous time transfer function as well as the discretized transfer functions using the various methods discussed. The continuous-time transfer function PID controller with a second order low-pass filter and two notches located at 950 and 1150 Hz. A sampling frequency of 5000 Hz was used.

The forward and backward Euler methods cannot track the controller correctly at around 300 Hz. Furthermore, they are unable to replicate the fast dynamical behavior of the two notches. As such, this is a clear example of why one should not use this methods. Tustin’s method shows the same dynamics as the continuous time transfer function, but due to the distortion in the frequency domain the notches are not located at the correct frequency. Tustin’s method combined with pre-warping is able to track the notches but loses magnitude when reaching the Nyquist frequency. In addition there is a phase offset. The Matched Pole Zero method is able to match the magnitude response very well, but has a large difference in phase which depending on the performance criteria might become problematic.

The Matlab code below was used to generate the plots.

clear all;
close all;
clc;

set(cstprefs.tbxprefs,'FrequencyUnits','Hz')

% Init variables and define continuous-time system
s = sym('s');
z = sym('z');
fs = 5e3;
h  = 1/fs;

Kp    = 1;
fd    = 100;
fi    = 10;
flp   = 400;
wlp   = 2*pi*flp;
blp   = 1;
 
fp{1} = 950;
fz{1} = 950;
bp{1} = 0.1;
bz{1} = 0.01;
wp{1} = 2*pi*fp{1};
wz{1} = 2*pi*fz{1};
 
fp{2} = 1150;
fz{2} = 1150;
bp{2} = 0.1;
bz{2} = 0.005;
wp{2} = 2*pi*fp{2};
wz{2} = 2*pi*fz{2};
 
Gc.sym = (wp{1}^2)/(wz{1}^2)*(s^2 + 2*bz{1}*wz{1}*s + wz{1}^2) / ...
                             (s^2 + 2*bp{1}*wp{1}*s + wp{1}^2) * ...
         (wp{2}^2)/(wz{2}^2)*(s^2 + 2*bz{2}*wz{2}*s + wz{2}^2) / ...
                             (s^2 + 2*bp{2}*wp{2}*s + wp{2}^2) * ...                     
         (Kp*(1/(2*pi*fd)*s + 1 + (2*pi*fi)/s)*(2*pi*flp)^2 / ...
         (s^2 + 2*blp*(2*pi*flp)*s + (2*pi*flp)^2));

[Gc.num,Gc.den] = numden(Gc.sym);
Gc.tf = tf(sym2poly(Gc.num),sym2poly(Gc.den));

% Forward Euler
Gd1.sym = subs(Gc.sym,s,(z-1)/h);
[Gd1.num,Gd1.den] = numden(Gd1.sym);
Gd1.tf = tf(sym2poly(Gd1.num),sym2poly(Gd1.den),h);

% Backward Euler
Gd2.sym = subs(Gc.sym,s,(z-1)/(h*z));
[Gd2.num,Gd2.den] = numden(Gd2.sym);
Gd2.tf = tf(sym2poly(Gd2.num),sym2poly(Gd2.den),h);

% Tustin
Gd3.sym = subs(Gc.sym,s,2*(z-1)/(h*(z+1)));
[Gd3.num,Gd3.den] = numden(Gd3.sym);
Gd3.tf = tf(sym2poly(Gd3.num),sym2poly(Gd3.den),h);

% Tustin with pre-warping
wlp   = 2/h*tan(wlp*h/2);
flp   = wlp/(2*pi);
wp{1} = 2/h*tan(wp{1}*h/2);
wz{1} = 2/h*tan(wz{1}*h/2);
wp{2} = 2/h*tan(wp{2}*h/2);
wz{2} = 2/h*tan(wz{2}*h/2);
      
Gc2.sym = (wp{1}^2)/(wz{1}^2)*(s^2 + 2*bz{1}*wz{1}*s + wz{1}^2) / ...
                              (s^2 + 2*bp{1}*wp{1}*s + wp{1}^2) * ...
          (wp{2}^2)/(wz{2}^2)*(s^2 + 2*bz{2}*wz{2}*s + wz{2}^2) / ...
                              (s^2 + 2*bp{2}*wp{2}*s + wp{2}^2) * ...                     
          (Kp*(1/(2*pi*fd)*s + 1 + (2*pi*fi)/s)*(2*pi*flp)^2 / ...
          (s^2 + 2*1*(2*pi*flp )*s + (2*pi*flp)^2));   

[Gc2.num,Gc2.den] = numden(Gc2.sym);
Gc2.tf = tf(sym2poly(Gc2.num),sym2poly(Gc2.den));

Gd4.sym = subs(Gc2.sym,s,2*(z-1)/(h*(z+1)));
[Gd4.num,Gd4.den] = numden(Gd4.sym);
Gd4.tf = tf(sym2poly(Gd4.num),sym2poly(Gd4.den),h);

% Matched Pole Zero
Gd5.tf = c2d(Gc.tf,h,'matched');

% Frequency spacing
f = logspace(0,log10(5000),3000);
wc = 2*pi*f;
wd = wc(wc < 2*pi*2500);

% Plot
figure;
hold all;
bode(Gc.tf,wc);
bode(Gd1.tf,wd);
bode(Gd2.tf,wd);
bode(Gd3.tf,wd);
bode(Gd4.tf,wd);
bode(Gd5.tf,wd);
xlim([1 3500]);
grid on;
legend({'Continuous-time','Forward Euler','Backward Euler', ...
    'Tustin','Tustin with pre-warping','Matched Pole Zero'});
Tags: , , , ,

Control System & PID controller

Control systems, as it is being taught on many universities in the world, regulates the behavior of a system using feedback. Systems may be mechanical, electrical, biological, economical or any other domain. Feedback as in many aspects of life is crucial. A target value r for an essential variable of the system, the process value y, is set. The target value is also referred as setpoint or reference value, whereas the process value is referred as measurement value. The control system calculates the difference between the setpoint and the process value, called the error e of the control system.

From this point we will focus mainly on the structure, design and implementation of the linear feedback controller. The figure below shows the structure of the control system with a linear feedback controller. We neglect any possible source of disturbance.

The linear feedback controller is often the combination a controller and multiple filters. The controller of the control system computes a control value u which is applied to the process and tries to regulate the error e to zero. A filter also often referred as compensator, is tries to improve the characteristics and performance.

The most famous feedback controller is the Proportional–Integral–Derivative (PID) controller is defined by (1). Based on the error it calculates a correction to be applied to the process based on proportional, integral, and derivative terms, respectively. Heuristically, the terms of a PID controller can be interpreted as corresponding to time: the proportional term depends on the present error, the integral term on the accumulation of past errors, and the derivative term is a prediction of future error, based on current rate of change.

(1)   \begin{equation*} H(s) = K_d \cdot s + K_p + K_i \frac{1}{s} \end{equation*}

Herein, K_d = K_p / (2\pi f_d) and K_i = K_p 2 \pi f_i for which f_d and f_i denote the frequency at which the derivative and integral action start to become active. In addition K_p denotes the proportional gain of the controller. An ideal derivative is not causal – its output depends on future inputs see (2) – therefore, the PID controller is often combined with either a first or a second order low-pass filter (3). This makes the PID controller causal, depending only on past and current inputs but not future inputs.

(2)   \begin{equation*} \frac{\mathrm{d}f(t)}{\mathrm{d}t} \triangleq \lim_{h \rightarrow 0} \frac{f(t + h) - f(t)}{h} \end{equation*}

(3)   \begin{equation*} H(s) = \left(K_d \cdot s + K_p + K_i \frac{1}{s}\right) \frac{\omega_{\text{lp}}^2}{s^2 + 2\beta\omega_{\text{lp}}s + \omega_{\text{lp}}^2} = \frac{\left(K_d s^2 + K_p s + K_i\right) \omega_{\text{lp}}^2}{s^3 + 2\beta\omega_{\text{lp}}s^2 + \omega_{\text{lp}}^2s} \end{equation*}

Note that the equation (3) results in a third order polynomial. With some algebra this can be reduced into a first order and second order polynomial. That is to say;

(4)   \begin{equation*} H(s) = \frac{K_i}{s} + \frac{\left(\omega^2_\text{lp} K_d - K_i\right)s + \left(\omega^2_\text{lp} K_p - 2\beta\omega_\text{lp} K_i\right)}{s^2 + 2 \beta \omega_\text{lp} s + \omega_\text{lp}^2}\end{equation*}

As such, the integral action is decomposed from the derivative and gain action along with the second order low-pass filter. With the integrator in parallel, it is possible to perform a saturation or anti-windup on the integrator signal.

Below you see an interactive bode diagram of a PID controller combined with a second order low-pass filter with some typical values. You clearly see the integrator action up to 10 [Hz], followed by the proportional gain action from 10 to 25 [Hz]. The derivative term is visible from 25 to 300 [Hz] and eventually the low-pass filter from 300 [Hz] and further.

Tags: , , ,