Lab02 PDF

Preview

Summary

Lab02...

Description

Control Engineering 3 - Laboratory Session 2 Objectives The goals of this lab are • Learn how to use MATLAB symbollic toolbox to take Laplace transform of a signal or function given in time-domain • Learn how to use MATLAB symbollic toolbox to take Inverse Laplace transform of a signal or function given in frequency domain. • Obtaining partial fraction expression for a transfer function using MATLAB • Defining transfer functions for Linear Time-Invariant(LTI) systems in MATLAB • Learning the ways to obtain poles and zeros of a transfer function in MATLAB • Learning to obtain step response for a given transfer function in MATLAB

Introduction to MATLAB Symbolic Toolbox MATLAB's Symbolic Math toolbox allows users to perform symbolic mathematical computations using MATLAB. The only basic requirement is to declare symbolic variables before they are used. For control systems analysis and design, symbolic math toolbox is particularly important because of the following: 1. Transfer functions and other expressions can be entered in symbolic form as we write them in the notebook. 2. Symbolic expressions can be manipulated algebraically and simplified. 3. It is straightforward to convert symbolic polynomials to the vectors of corresponding power-term coefficients. Type help symbolic on the MATLAB command prompt to see all the functionalities available with Symbolic Math toolbox. In this section, we will learn commands of Symbolic Math toolbox particularly useful for control engineering.

Example

syms f(x) f(x) = x^4-2*x^3+6*x^2-2*x+10 f(x) =

1

will give you the function in terms of variable x. For example, to evaluate this function at f(-5) ans =

Example Find the intersection between lines

and

using solve. Equate the lines using the "==" operator.

syms y1 y2 y1 = x+3; y2 = 3*x; solve(y1 == y2) ans =

The Symbolic Math Toolbox supports evaluation of mathematical functions by substituting for any part of an expression using subs. You can substitute numeric values, other symbolic variables or expressions, vectors, or matrices. The Symbolic Math Toolbox supports the solving of equations and systems of equations using solve. It supports solving multivariate equations, solving inequalities and solving with assumptions. Solutions can be found symbolically or numerically with high precision by using variable-precision arithmetic.

Example Make substitutions with your symbolic variables. Substitute

into

syms x xo subs(x^2+1,x,xo-1) ans =

Here the second argument of subs command is the old variable and the last one is the new variable. Note that subs(s,old,new) returns a copy of s, replacing all occurrences of old with new, and then evaluates s.

2

Example For example, evaluate

by substituting

,

,

syms a b c subs(cos(a) + sin(b) - exp(2*c), [a b c], [pi/2 pi/4 -1]) ans =

Example Find the zeros(roots) of the polynomial equation

.

syms x solve(9*x^2 - 1 == 0) ans =

Question Solve the general quadratic equation for . , ,

and use subs to evaluate that solution

We can also multiply polynomials and simplify the results easily in MATLAb. For instance

Example Perform the following polynomial multiplication and simplify the results (by using simplify command):

simplifies to

3

simplify((x - 1)*(x + 1)*(x^2 + x + 1)*(x^2 + 1)*(x^2 - x + 1)*(x^4 - x^2 + 1)) ans =

Question Find the roots(zeros) of the following polynomial given in multiplication form

syms x solve((x^2 + 3*x + 2)*(x^2 - 9)==0) ans =

We can also factor or expand multivariate polynomials in MATLAB Symbolic toolbox.

Example syms x y factor(y^6-x^6) ans =

Laplace Transform One of the most useful properties for people in control is dealing with Lapalce transform under Symbolic toolbox. laplace(f) returns the Laplace Transform of f. Recall that

4

where

is a complex variable. Thus knowing

function

, that is called the Laplace Transform of

Conversely, when

and

.

given, the inverse Laplace Transform of

allows us to get the function of

and the integral given above exists, we can find a

where

,

is a unit step function and multiplication

yields a time function that is zero for

.

However, taking the Laplace transforma or Inverse Laplace Transform of a general function by hand is not that easy unless some simple functions are considered.Therefore, one can use MATLAB Symbolic Toolbox to dela with these transformations. Let us see how it works

Example Compute the Laplace transform

.

syms t laplace(heaviside(t)) ans =

Example Compute the Laplace transform

.

laplace(heaviside(t-5)) ans =

5

Here the heaviside function works as a unit step function. ezplot(heaviside(t-5), [-5, 10])

Example Consider a time domain signal charectrised by

Its Laplace transform can be obtained as syms t s f = (-1.25+3.5*t*exp(-2*t)+1.25*exp(-2*t)) f =

6

F=laplace(f,t,s) F =

One can obtain a more tidy expression by using simlify command. Hence, simplify(F) ans =

Inverse Laplace Transform The command one uses now is ilaplace. One also needs to define the symbols t and s. Lets calculate the inverse of the previous function syms t s F=(s-5)/(s*(s+2)^2); ilaplace(F) ans =

Partial Fractions in MATLAB To find the inverse Laplace transform of a complicated function, we can convert the function to a sum of simpler terms for which we know the Laplace transform of each term. The result is called the partial-fraction expansion. If , where the order of is less than or equal the order of . The residue function converts a quotient of polynomials to pole-residue representation, and back again to polynomials. [r,p,k]=residue(B,A) 7

finds the residues, poles and direct term of a partial fraction expansion of the ratio of two polynomials . The quotient of polynomials and is represented as

where is a column vector of residues, is a column vector of pole locations (the roots of the denominator), and k is a row vector of direct terms (if the degree of the numerator is less than the degree of the denominator, then k will be an empty vector). The residue command requires two input vectors: one holding the coefficients of the numerator and one holding the coefficients of the denominator. The coefficients must be given in descending powers of s. Consider for instance the following function:

To find the partial fraction decomposition of the denominator as vectors. The numerator,

, we must first enter the polynomials at the numerator and at , can be entered in MATLAB as

N = [5,-1] % coefficients of the numerator in decreasing order N = 1×2 5

-1

while the denominator is represented by the vector D = [1, 0, -3, -2] % coefficients of the denominator in decreasing order D = 1×4 1

0

-3

-2

Note that the term is missing, thus we need to enter a zero as the second entry of the vector. The residue command gives the following output: [r, p, k] = residue(N,D) r = 3×1 1.0000 -1.0000 2.0000 p = 3×1 2.0000 -1.0000 -1.0000 k = []

Thus the partial fraction decomposition of

is 8

Please note that The last term in the partial fraction decomposition follows from the fact that the pole −1.0000 is repeated.

Solving ODEs using the Laplace and Inverse Laplace Transform Example A mass-spring system is modeled by the equation

where Our task is to find the output trajectory First, Use MATLAB to find the Laplace transform of the forcing term unit step function is denoted in MATLAb by heaviside.

We can write the mathematical expression for

and plot the forcing term. Note: The

as follows:

To see how this forcing term varies in time and to obtain its Laplace transform, one can use the following commands: syms s t g = 8*t + heaviside(t-5)*(40-8*t); fplot(g,[0,10])

9

G= laplace(g) G =

G = collect(G) G =

Here fplot(FUN,LIMS) plots the function FUN between the x-axis limits specified by LIMS = [XMIN XMAX]. Applying the Laplace transform to both sides of the Differential Equation and substituting the Initial Conditions gives

10

where

is the Laplace transform of the solution

. Let

We need to use residue command to find the partial fraction decomposition of F ans use this partial decomposition to find (by hand) the inverse Laplace transform of . num=[8]; % coefficients of the numerator in decreasing order den=[1,0,4,0,0]; % coefficients of the denominator in decreasing order [r,p,k]=residue(num,den) r = 4×1 complex 0.0000 + 0.5000i 0.0000 - 0.5000i 0.0000 + 0.0000i 2.0000 + 0.0000i p = 4×1 complex 0.0000 + 2.0000i 0.0000 - 2.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i k = []

Thus

and combining the complex terms we obtain

Using a Table of Laplace transforms, we find that the inverse Laplace transform of

and the inverse laplace of

is

where the response can be plotted as follows: y = 2*t-sin(2*t)-heaviside(t-5)*(2*(t-5)-sin(2*(t-5))) y =

11

is

fplot(y, [0 10]) Let us see both graphs on the same frame: hold on fplot(g,[0,10]) xlabel('time'); ylabel('$g(t)$, $y(t)$', 'Interpreter','latex'); hold off

Here we used xlabel and ylabel commands to put labels for the x-axis and y-axis, respectively. The command hold on is used to retains plots in the current axes so that new plots added to the axes do not delete existing plots. Similarly, hold off sets the hold state to off so that new plots added to the axes clear existing plots We can see that the graph of the solution makes perfect sense in the mass-spring context. The linearly increasing forcing term for the first 5 seconds causes the equilibrium to shift linearly. Once the forcing term

12

stabilizes into a constant value so does the equilibrium and, since the system is undamped, the mass keeps oscillating around the new equilibrium with constant amplitude.

Question Consider the transfer function

1. Use the MATLAB function residue to find the residues and poles of

. Use the output to find the

partial fraction decomposition of . 2. Use the Table of Laplace transforms and your answer to (1) to find the inverse Laplace transform of . 3. Enter the function Y in MATLAB and confirm your result from part (2) with the ilaplace command (do not forget to declare s as a symbolic variable).

Question Confirm both by hand and computer that the Laplace transform of

is

Question A mass is attached to a spring with constant and damping constant . The mass is initially in equilibrium position and it is released with an initial velocity of 2 units. At the instant , the mass is struck with a hammer and at t it is struck again. The motion of the mass satisfies the following differenetial equation

where the initial conditions are

and

. The delta (dirac) function is entered in MATLAB with

the command dirac. • Use MATLAB to find the Laplace transform of the right hand side. • Apply (by hand) the Laplace transform to both sides of the equation and determine transform of the solution

. Use ilaplace to find the inverse Laplace of

13

, the Laplace

• Plot the solution in the interval and in relation to the forcing term.

, and comment on the graph in the mass-spring context

Transfer functions and time responses in MATLAB A Transfer Function is the ratio of the output of a system to the input of a system, in the Laplace domain considering its initial conditions and equilibrium point to be zero. You can create transfer function (TF) models by specifying numerator and denominator coefficients. For example, to define

one can use the following lines: num = [1 0]; den = [1 3 4]; G = tf(num,den) G = s ------------s^2 + 3 s + 4 Continuous-time transfer function.

A useful trick is to create the Laplace variable, s. That way, you can specify polynomials using s as the polynomial variable. For example s = tf('s'); G = s / (s^2 + 3*s + 4) G = s ------------s^2 + 3 s + 4 Continuous-time transfer function.

One can also Create Transfer Function Model Using Zeros, Poles, and Gain. This example shows how to create single-input, single-output (SISO) transfer functions in factored form using zpk,

Example Create the factored transfer function

14

Z = [0]; P = [-1-1i -1+1i -2]; K = 5; G = zpk(Z,P,K); Z and P are the zeros and poles (the roots of the numerator and denominator, respectively). K is the gain of the factored form.

System Interconnections by MATLAB command line tools The Matlab control system toolbox has various functions to connect LTI models. In the following section we introduce the three functions series, parallel and feedback. These can be very useful for creating larger systems based on smaller subsystems, for example a plant with its controller.

Series Interconnection This command connects two LTI models in series (cascade). the syntax of the command: sys = series(sys1,sys2) series connects two LTI models in series. This function accepts any type of LTI model. This command is equivalent to the direct multiplication sys = sys2 * sys1.

Example Assume

and

. to connect sys1 and sys2 in series we need to write the following

commands: s=tf('s'); sys1 = 1/(s+1); sys2 = (s+1)/(s*(s+3));

15

sys = series(sys1,sys2) sys = s + 1 ----------------s^3 + 4 s^2 + 3 s Continuous-time transfer function.

Please note that the the series command does not perform any cancellations between zeros and poles of the resultant transfer function. Therefore we need to perform these cancellations by using a separate command minreal as follows: minreal(sys) ans = 1 --------s^2 + 3 s Continuous-time transfer function.

Parallel Interconnection Parallel interconnection is performed by the command parallel. and the syntax of the command is sys = parallel(sys1,sys2) This command is equivalent to the direct addition sys = sys1 + sys2.

Example Assume

and

. To connect sys1 and sys2 in parallel we need to write the following

commands:

16

s=tf('s'); sys1 = 1/(s+1); sys2 = (s+1)/(s*(s+3)); sys = parallel(sys1,sys2) sys = 2 s^2 + 5 s + 1 ----------------s^3 + 4 s^2 + 3 s Continuous-time transfer function.

Feedback interconnection This command applies a feedback interconnection on two LTI models. sys = feedback(sys1,sys2) returns an LTI model sys for the negative feedback interconnection. The closed-loop model sys has u as input vector and y as output vector. Note that to apply positive feedback, use the syntax sys = feedback(sys1,sys2,+1)

Example Assume

and

. To connect sys1 and sys2 in feedback structure with a newgative

feedback, we need to write the following commands: s=tf('s'); sys1 = 1/(s+1); sys2 = (s+1)/(s*(s+3));

17

sys = feedback(sys1,sys2) sys = s^2 + 3 s --------------------s^3 + 4 s^2 + 4 s + 1 Continuous-time transfer function.

sys = minreal(sys) % Again we need to use minreal to perform any pole-zero cancellations that sys = s^2 + 3 s --------------------s^3 + 4 s^2 + 4 s + 1 Continuous-time transfer function.

Questions For the following multi-loop feedback system, obtain closed loop transfer function

first

by hand then by MATLAB command based tools and verfiy your answer. Note: Do not forget to use minreal command for simplifications.

where ,

,

,

,

and

Step response step calculates the step response of a dynamic system. step(sys,Tfinal) simulates the step response from t = 0 to the final time t = Tfinal. Express Tfinal in the system time units, specified in the TimeUnit property of sys. 18

Example To obtain the step response of a unity feedback system

with

, one can use the following set of commands

s = tf('s'); G = exp(-s) * (0.8*s^2+s+2)/(s^2+s); Gcl =feedback(G,1); step(Gcl)

19

Question Consider the feedback system depicted in the figure below

• Compute the closed-loop transfer function using the ‘series’ and ‘feedback’ functions • Obtain the closed-loop system unit step response with the ‘step’ function and verify that final value of the output is 2/5.

Question A satellite single-axis altitude control system can be represented by the block diagram in the figure given below.

20

The variables k, a and b are controller parameters, and J is the spacecraft moment of inertia. Suppose the nominal moment of inertia is , and the controller parameters are , , and . •

Develop an m-file script to compute the closed-loop transfer function

• Compute and plot the step response to a step input. • The exact moment of inertia is generally unknown and may change slowly with time. Compare the step response performance of the spacecraft when Jis reduced by 20% and 50%. Discuss your results.

Finally the command stepinfo computes the step response characteristics. the syntax of the command is as follows: S = stepinfo(Y,T,YFINAL) takes step response data (T,Y) and a steady-state value YFINAL and returns a structure S containing the following performance indicators: * RiseTime: rise time * SettlingTime: settling time * SettlingMin: min value of Y once the response has risen * SettlingMax: max value of Y once the response has risen * Overshoot: percentage overshoot (relative to YFINAL) * Undershoot: percentage undershoot * Peak: peak absolute value of Y * PeakTime: time at which this peak is reached.

For example to obtain the rise-time of the step response of Gcl, one can use the following set of commands S = stepinfo(Gcl); S.RiseTime ans = 0.2464

21

and for the settling time: S.SettlingTime ans = 16.0130

22...