Archive for category Chapter 3. Closed loop TF

Slides session 10

The last bunch of slides:

ppt Session10_Week12

Leave a comment

Semilog paper

The slide below is an example of a semilog grid, that can be used to draw Bode plots by hand.

ppt Semilog

Leave a comment

Slides session 9

Exercise on time delay in the loop:

: ppt Session9_Week11

Matab code:

% Bode exercise destillation
% pkg load control; % Octave only
close all

% Process without time delay
Tpu = tf(1,[10 1]);
% Control parameters
Kc = 10;
TI = 2;
TD = 0;
taud = 0; % filter for D action
% Controller
Tc = Kc*(1+tf(1,[TI 0])+tf([TD 0],[taud 1]))
% Open loop without time delay
sysol = Tpu*Tc;
bode(sysol)

% Time delay
L = 1;
% Creating an approximate TF for the time delay since Octave doesn't recognize e^s
[num,den] = pade(L,15);
systd=tf(num,den);
% Open loop with time delay
sysoltd = sysol*systd;
% Closed loop
syscltd = feedback(sysoltd,1)

% Bode plot of the open loop with time delay
figure
bode(sysoltd,[0.01:0.1:10],'r')
% Step response of the closed loop
figure
step(syscltd)

This Matlab code uses the Pade command, which is not as such available in Octave, so I’m not sure it works the same in Matlab. So just for certainty, here’s pade.m

# Copyright (C) 2013 Thomas Vasileiou
##
## This file is part of LTI Syncope.
##
## LTI Syncope is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## LTI Syncope is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {} pade (@var{T}, @var{N})
## @deftypefnx {Function File} {[@var{num}, @var{den}] =} pade (@var{T}, @var{N})
## @deftypefnx {Function File} {@var{retsys} =} pade (@var{sys}, @var{N})
## @deftypefnx {Function File} {@var{retsys} =} pade (@var{sys}, @var{Nu}, @var{Ny})
## @deftypefnx {Function File} {@var{retsys} =} pade (@var{sys}, @var{Nu}, @var{Ny}, @var{Nint})
## Time delay exp(-T*s) approximation by rational function using
## the Pade expansion.
##
## If the time of the delay @var{T} and the order of the approximation
## @var{N} are provided, @code{pade} returns the nominator @var{num} and
## the demonimator @var{den} of the expansion. If the function is called
## without output arguments, a comparison plot of the step response and 
## the phase between the pure delay and the approximation is shown
## instead.
## 
## Direct substitution of the delays of an @acronym{LTI} system is possible
## by passing the system directly to the function. For discrete time
## models, this function redirects to @code{absorddelay (@var{sys})}.
## 
##
## @strong{Inputs}
## @table @var
## @item sys
## Input system with delays.
## @item T
## Time delay.
## @item N
## Order of the Pade approximation. For models with multiple delays,
## using a single order @var{N} will result in approximating all
## delays by the same order.
## @item Nu, Ny
## For models with multiple channels and delays, different orders of
## approximation may be achieved. Use the vector @var{Nu} to provide the
## order of each input channel and the vector @var{Ny} to provide the order
## for each output channel. If a single scalar is provided for
## @acronym{MIMO} systems, the same order is used for all input or output
## channels, respectivelly. If any of provided order is an empty matrix,
## the delays related to corresponding channel are not approximated. The value
## @code{inf} can be used in the vectors to avoid approximation on specific
## channels, whereas the value 0 just removes the delay (0th-order approximation).
## @item Nint 
## The matrix @var{Nint} is used to determine the orders of approximation
## of the io or the internal delays of the system. If a single scalar is 
## provided for @acronym{MIMO} systems or for ss models with multiple intenal
## delays, the same order is used for all delayed channels. If this argument
## is ommited or an empty matrix is provided, io or internal delays are not
## approximated. The value @code{inf} can be used in the matrix to avoid
## approximation on specific channels, whereas the value 0 just removes the
## delay (0th-order approximation).
## @end table
##
## @strong{Outputs}
## @table @var
## @item num, den
## Numerator and denominator of the expasion of exp(-T*s).
## @item retsys
## Resulting system.
## @end table
##
## @strong{Example}
## @example
## @group
## octave:1> [num, den] = pade (0.2, 2)
##
## num =
##
## 0.040000 -1.200000 12.000000
##
## den =
##
## 0.040000 1.200000 12.000000
## @end group
## @end example
## @example
## @group
## octave:2> mm1 = tf (1 ,[1 1]);
## octave:3> mm1.indelay = 0.2;
## octave:4> mm1.outdelay = 0.2;
## octave:5> mm1.iodelay = 0.3
##
## Transfer function 'mm1' from input 'u1' to output ...
## 
## 1 
## y1: e^(-0.7 s) -----
## s + 1
## 
## Continuous-time model.
## @end group
## @end example
## @example
## @group
## octave:6> pade (mm1, 1)
##
## Transfer function 'ans' from input 'u1' to output ...
##
## -0.012 s^3 + 0.32 s^2 - 2.8 s + 8 
## y1: ---------------------------------------------
## 0.012 s^4 + 0.332 s^3 + 3.12 s^2 + 10.8 s + 8
##
## Continuous-time model.
## @end group
## @end example
## 
## @example
## @group
## octave:6> pade (mm1, inf, 2, 0)
## 
## Transfer function 'ans' from input 'u1' to output ...
## 
## 0.04 s^2 - 1.2 s + 12 
## y1: e^(-0.2 s) ---------------------------------
## 0.04 s^3 + 1.24 s^2 + 13.2 s + 12
## 
## Continuous-time model.
## @end group
## @end example
##
## @seealso{absorbdelay}
## @end deftypefn

## Author: Thomas Vasileiou <thomas-v@wildmail.com>
## Created: January 2013
## Version: 0.1

function [num, den] = pade (sys, nu, ny, nint)
 
 ## arguments checking
 lti_flag = isa (sys, "lti");
 
 ## if not an lti
 if (nargin == 2 && ! lti_flag)
 ## no sys, calculate coeffs
 if (! is_real_scalar (sys) || sys < 0)
 error ("pade: Parameter T must hold a real positive scalar value.")
 endif 
 
 if (! is_real_scalar (nu) || nu < 0)
 error ("pade: The order N must be a real non-negative scalar value.")
 endif
 nu = round (nu);
 
 ## if no output argument, do the ploting
 if (nargout == 0)
 plot_pade (sys, nu);
 else
 [num, den] = padecoeff (sys, nu);
 endif
 return;
 endif
 
 if (! lti_flag)
 error ("pade: when called with more than two arguments, the first one must be an LTI system.")
 endif
 
 ## general parameters
 [p, m] = size (sys);
 if (isa (sys, "ss"))
 nint_ss = length (get (sys, "internaldelay"));
 else
 nint_ss = [];
 endif
 
 ## no delay? return
 if (! hasdelay (sys))
 num = sys;
 return;
 endif
 
 switch (nargin)
 case 2
 if (! is_real_scalar (nu))
 error ("pade: When providing a single order N for all delays, it must be a scalar value.")
 endif
 [nu, ny, nint] = __check_pade_orders__ (nu, nu, nu, p, m, nint_ss);
 
 case 3
 [nu, ny, nint] = __check_pade_orders__ (nu, ny, [], p, m, nint_ss);
 
 case 4
 [nu, ny, nint] = __check_pade_orders__ (nu, ny, nint, p, m, nint_ss);
 
 otherwise
 print_usage (); 
 endswitch
 
 num = __pade__ (sys, nu, ny, nint);
 den = [];

endfunction


function [delay_num, delay_den] = padecoeff (t, n)

 if (t > 0)
 ## calculate pade coefficients as:
 ##
 ## (2*N - j)!
 ## For order N: a_j = ----------------- 
 ## j! * (N - j)! 
 aj = zeros (1, n+1);
 for ij = 1 : n+1
 aj(ij) = factorial (2*n-ij+1)/factorial (n-ij+1)/factorial (ij-1);
 endfor
 

 ## add the time
 delay_num = zeros (1, n+1);
 delay_den = zeros (1, n+1);
 for ij = 1 : n+1
 delay_den(ij) = aj(n+2-ij)*t^(n+1-ij);
 delay_num(ij) = (-1)^(n+1-ij)*delay_den(ij);
 endfor
 else
 delay_num = 1;
 delay_den = 1;
 endif
 
endfunction


function plot_pade (t, n)

 ## get coeffs & make system
 [delay_num, delay_den] = padecoeff (t, n);
 delay_sys = tf (delay_num, delay_den);
 
 ## step and phase response 
 [yd, td] = step (delay_sys, linspace (0, 2*t, 100));
 [~, phd, wd] = bode (delay_sys);
 
 ## make starting from 0 deg
 if (phd(1) > 0)
 phd = phd - 360;
 endif

 #figure;
 subplot (2, 1, 1);
 plot ([0, t, t, 2*t], [0, 0, 1, 1], '--r');
 hold on;
 plot (td, yd);
 hold off;
 title (['Step Response of ', num2str(n), '-order Pade approximation']);
 ylabel ('Amplitude');
 xlabel ('Time [s]');
 legend ('pure delay', 'Pade approximation', 'Location', 'SouthEast');
 yaxis_min = 1.4*min (yd);
 yaxis_max = 1.2*max (yd);
 if (yaxis_min >= yaxis_max)
 yaxis_min = 0;
 yaxis_max = 1;
 endif
 axis ([0, 2*t, yaxis_min, yaxis_max]);

 subplot (2, 1, 2);
 semilogx (wd, -wd*180/pi*t, '--r');
 hold on;
 semilogx (wd, phd);
 hold off;
 title ('Phase response comparison');
 ylabel ('Phase [deg]');
 xlabel ('Frequency [rad/s]');
 axis ([wd(1), wd(end), min(phd) - 500, 0]);

endfunction


function [nu, ny, nint] = __check_pade_orders__ (nu, ny, nint, p, m, nint_ss)

 ## chech input orders 
 if (! isempty (nu))
 if (! is_real_vector (nu) || any (nu + eps < 0))
 error ("pade: parameter 'Nu' requires an integer positive vector or scalar.");
 endif
 nu = ceil (nu);
 if (is_real_scalar (nu))
 nu = nu*ones (m, 1);
 elseif (length (nu) != m)
 error ("pade: parameter 'Nu' requires a vector with as many elements as the inputs of the system.");
 else
 nu = reshape (nu, [], 1); 
 endif 
 endif
 
 ## check output orders 
 if (! isempty (ny))
 if (! is_real_vector (ny) || any (ny + eps < 0))
 error ("pade: parameter 'Ny' requires an integer positive vector or scalar.");
 endif
 ny = ceil (ny);
 if (is_real_scalar (ny))
 ny = ny*ones (p, 1);
 elseif (length (ny) != p)
 error ("pade: parameter 'Ny' requires a vector with as many elements as the outputs of the system.");
 else
 ny = reshape (ny, [], 1);
 endif 
 endif
 
 ## check iodelay orders 
 if (! isempty (nint))
 if (isempty (nint_ss))
 ## not an ss model, check io matrix 
 if (! is_real_matrix (nint) || nnz (nint + eps < 0))
 error ("pade: parameter 'Nint' requires an integer positive matrix or scalar.");
 endif
 nint = ceil (nint);
 if (is_real_scalar (nint))
 nint = nint*ones (p, m);
 elseif (! isequal (size (nint), [p, m]))
 error ("pade: parameter 'Nint' requires a matrix with as many columns as the inputs and as many rows as the outputs of the system.");
 endif 
 else
 ## ss model
 if (nint_ss > 0)
 ## has internal structure
 if (! is_real_vector (nint) || any (nint + eps < 0))
 error ("pade: parameter 'Nint' requires an integer positive vector or scalar.");
 endif
 nint = ceil (nint);
 if (is_real_scalar (nint))
 nint = nint*ones (nint_ss, 1);
 elseif (length (nint) != nint_ss)
 error ("pade: parameter 'Nint' requires a vector with as many elements as the internal delays of the system.");
 else
 nint = reshape (nint, [], 1);
 endif 
 endif
 endif
 endif
 
endfunction


## TODO: add more tests for multiple channels ... (and maybe state space)
%!shared num, den, expnum, expden
%! expden = [1, 1.669, 1.692, 0.7692];
%! expnum = [-1, 0.7692];
%! sys = tf (1, [1, 0.9, 1], "outdelay", 2.6);
%! sys = pade (sys, 1);
%! [numt, dent] = tfdata (sys, "vector");
%! den = dent / dent(1);
%! num = numt / dent(1);
%!assert (num, expnum, 1e-3);
%!assert (den, expden, 1e-3);
%! expden = [1, 5.515, 14.03, 19.43, 15.02, 6.827];
%! expnum = [-1, 4.615, -8.876, 6.827];
%! sys = tf (1, [1, 0.9, 1], "outdelay", 2.6);
%! sys = pade (sys, 3);
%! [numt, dent] = tfdata (sys, "vector");
%! den = dent / dent(1);
%! num = numt / dent(1);
%!assert (num, expnum, 2e-3);
%!assert (den, expden, 2e-3);

%!error (pade (-1, 1));
%!error (pade (1, -1));

Leave a comment

Slides session 8

(Session 7 – Bode plots: no slides, see Appendix)

Gain and phase margin/time delays (session by P. Spaepen):

ppt Session8_Week10

Leave a comment

Slides session 6

Closed loop transfer functions:

ppt Session6_Week7

and here’s some Matlab code that goes along with the exercises in the slides:

% Closed loop TF
% pkg load control;
close all

% Oef 1 en 2
Tpu = tf(8,[1 8 4]);
Kc = 1.5;
syscl = feedback(Kc*Tpu,1);
step(syscl);
hold on
Kc2 = 7.5;
syscl2 = feedback(Kc2*Tpu,1);
PO = 100*exp(-0.5*pi/sqrt(1-0.5^2))
Tpeak = 2*pi/8/sqrt(1-0.5^2)
step(syscl2)

% Oef 3
Tpu3 = tf(2,[1 1 1]);
Kc3 = -3/8;
syscl3 = feedback(Kc3*Tpu3,1);
step(syscl3);

%Oef4
Tpu4 = tf(1,[1 4 3]);
Ts4 = tf(1,[1 5]);
syscl4 = feedback(192*Tpu4,Ts4)
step(syscl4,10)

And here’s a nice video with some nice exercises that looks at the closed loop TF from a slightly different angle (forward path/(1+return path)).
Watch out, in the first exercise, the question is to find the closed loop TF from reference to u (whereas up so far, we only treated the r to z case).

Leave a comment

Slides session 5

Slides from session 5 on transfer functions.

 Session5_Week6

Leave a comment

Follow

Get every new post delivered to your Inbox.

Join 241 other followers