# Chapter-Title by sdaferv

VIEWS: 17 PAGES: 12

• pg 1
```									MATLAB for Systems Analysis
We know how to model systems, and MATLAB. Let’s put these together. We will consider the use of MATLAB for obtaining Response in time domain Response in frequency domain For the former you will be given help to allow you to write code to Find the theoretical response to a step input Find the simulated response to a step input For the latter you will be given help to allow you to develop code to Plot Nyquist and Bode diagrams of the loop gain (strictly -loop gain) Plot Frequency response of whole, or closed loop system The assignment associated with this work will involve various systems. These systems will be represented using transfer functions, which are in turn represented using polynomials.
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.1

Example System
I controller plant 1 1+s15 9 (1+s15)(1+s8) 1+s10 O

(Minus) Loop transfer function is 1  s15 1 9 9 9 *   1  s10 (1  s15)(1  s8) (1  s10)(1  s8) 80s 2  18s  1 Closed Loop transfer function 9 9 9 80s 2  18s  1   9 80s 2  18s  1  9 80s 2  18s  10 1 80s 2  18s  1
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.2

In MATLAB, represent transfer function with polynomials: c = tf (9*[15 1], [10 1]); p = tf (1, conv([15 1],[8 1]) ); For some frequency response need loop transfer function: series (c,p); % multiplies two transfer functions in series This does not do cancelling, so simplify by: minreal (series(c,p)); For closed loop analysis, want closed loop transfer function: minreal (feedback (series (c,p), 1) ); % feedback (t1,t2) is t1/(1+t1*t2). MATLAB code to generate transfer functions of other components Trans Func

5s

2

3s 2  2s  1 MATLAB tf([5 0],1); tf(2, [3 2 1]); tf([4 5], [1 0]); tf((1,4,5),[1 0]) Note, for those without the control toolbox, I have provide versions of the above MATLAB functions: these are on a separate sheet.
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.3

5 4 s

1 4   5s s

Theoretical Response
The Laplace transform lectures mean you should be able to find the K theoretical step response for a system with transfer function . 2  bs  c as The following is part of a MATLAB function to return such a step response function ans = sotheory(k,a,b,c,t); % returns theoretical response of system with tf k/(as^2+bs+c) % at all the values of time in the vector t. if (a == 0) & (b ~= 0) & (c ~= 0) ans = (k/c)* (1 - exp(-t*c/b)); elseif (a ~= 0) & (b == 0) & (c ~= 0) w = sqrt(c/a); ans = k*(1-cos(w*t))/c; end NB to evaluate exp(dt)sin(wt) do exp(d*t).*sin(w*t) Note the '.*' operator
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.4

Step Response of Such a System
The command step returns the simulated step response and plots it. It needs the numerator and denominator polynomials, with optional time. Can call with output arguments, to get data but not plot it. The following code sets values for t, uses sotheory to get the theoretical response, plots that response, and then calls step to add the simulated response. It assumes values for variables k, a, b and c exist. t = [0:0.1:50]; % set range of values of t plot (t, sotheory(k, a, b, c, t) ); % calc and plot theoretical resp hold on; % hold data for the graph step (k,[a b c],t, ‘*’); % calc and plot step response Alternative method [op, dummy, t] = step(k, [a b c]); % use step to calc values and t plot (t, op, t, sotheory (k, a, b, c, t), ‘*’); % plot step & calc & plot theory % note dummy variable (it is the 'state variable', which we can ignore here) % note also, plot one graph in * so can see agreement of two methods
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.5

O 1  sTe As an intro, consider simulating system with transfer function  I 1  sTa x 1 O Could do by:  and  1  sTe Here x is ‘state’ variable. I 1  sTa x Conceivably do by:

s I 1 Ta 1 X s

Te

O

X(n) - X(n - 1) But, this involves differentiating, done by t n - t n -1 X values have errors, and tn-tn-1 is small, so errors amplified. So avoid.
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.6

Te I 1 Ta 1 X s

O

Uses differential of X we have already got. Thus, algorithm to find O at a given time is: Sx = (I – X)/Ta; % find input to integrator X = Integrate (Sx); % integrate this to get X O = X + Sx * Te; % compute O Can extend concept to system with transfer function O z n s n z n -1s n -1 .. z1s z 0 Note, pn <> 0.  I p s n p s n -1 .. p s p n n -1 1 0
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.7

Zn Zn-1 Z1 Z0

S

O

I

1 Pn

1 s

1 s

1 s

x
Po P1 Pn-2 Pn-1

S

2/CY/I6 Modelling project

MATLAB System Analysis

RJM p 5.8

We again define state variable x, saying x 1 O  and  z n s n z n -1s n -1 .. z1s z 0 I p s n p s n -1 .. p s p x
n n -1

From x/I we get

And O/x gives Note, we find nth differential of x and by integrating get n-1th differential, etc. To find O and x we need these differentials, but we already have them. Note can find values of O, x and its diffs at any time with vectors: num = numerator poly, den = den poly, and x is [s nx sn-1x …sx x] % if a=[a1 a2 a3] and b=[b1 b2 b3], then dot(a,b) = a1*b1 + a2*b2 + a3*b3 x(1) = (Input - dot (den(2:n), x(2:n)) / den(1); % find nth diff of x Integrate x(2:n) % integrate to find n-1th, etc. Output = dot(x,num) % compute output
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.9

1 0 n x  1 I - p s n -1x .. p sx p x s n -1 1 0 pn O  z n s n x z n -1s n -1x .. z1sx z 0 x





Problem - how to integrate
Recall integration is area under a curve
Previous Input 1 s Current Input

Input

Output Previous Time Current Time

Assume we know integral up to previous time. Then Euler's (poor) rectangular method is Area = Area + Area of Rectangle due to current input = Area + (CurrentInput) * (CurrentTime - Previous Time) Note Runge-Kutta is better integration method. Hence we have the following MATLAB function to find the step response.
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.10

function resp = RJMStep (num, den, t) % Euler Simulation of system with transfer function num/den, at times in t % store output response in resp. Note den(1) must not be zero. % Dr Richard Mitchell, 19/9/00 x = zeros (size(den)); % x is array of state variables: set to 0 n = length(den); % find length of den array while length(num) < n % add leading zeros to num, so same length as den num = [0 num]; end; resp = [0]; % initialise response array to 0 for ct = 2:length(t), % for all time x(1) = (1 - dot(den(2:n),x(2:n)))/den(1); % calc input to first integrator x(2:n) = x(2:n) + (t(ct)-t(ct-1))*x(1:n-1); % integrate all state vars resp = [resp dot(num,x)]; % add output to resp array end; Note to be reasonably accurate, t(ct) - t(ct-1) must be sufficiently small.
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.11

Example:

5 s 2  2s  5
1.4 1.2 1 0.8 0.6 0.4 0.2 0

if t = [0:0.1:10], t=[0:0.02:10] and SOTheory

0

1

2

3

4

5

6

7

8

9

10

Not bad if Tstep = 0.02, but more calcs. Need better integration method
2/CY/I6 Modelling project MATLAB System Analysis RJM p 5.12

```
To top