ChEE Department    IIT

Matlab Tutorial for Process Control


 


Matlab is an abbreviation for MATrix LABoratory. Matlab is a high-level programming environment that processes arrays and matrices and provides a powerful graphical environment. A high-level programming environment allows the users to program without worrying about declaring variables, allocating memory, using pointers, and compiling code and other routine tasks, which are associated with languages such as FORTRAN and C. Matlab also incorporates many built in functions that can perform a variety of complex mathematical routines, from finding eigenvalues to solving differential equations. The purpose of this tutorial is to familiarize the user to working in Matlab. Additional information can be found by using Matlab’s extensive help library. Detailed information about Matlab commands and functions may be obtained in three different ways: (i) by calling help files from the Help button in the Matlab command window or by typing help at the Matlab prompt (»); (ii) by reading the html files via a browser (you may need to activate your browser first); (iii) by calling the Online Manuals from the Help Desk Main Menu that is displayed by the browser. The following subjects will be covered in this tutorial with exercises following certain sections.

Sections

  1. Starting Matlab
  2. The Command Line and Basic Commands
  3. Script Files and Functions
  4. Using Vectors and Matrices
  5. Plotting and Graphics
  6. Programming in Matlab
  7. Solving Ordinary Differential Equations
  8. Process Control Functions
  9. Closed-loop Model Development Example: Mixing Process Model

  10.  
This tutorial is more effective when the reader is interacting with Matlab. Please go to the PC lab in Perlstein Hall and work through it!

Written by Jeff DeCicco August 31, 1994
Revised by Jeff DeCicco January 30, 1995
Revised by Ali Cinar September 9,1998
 

1) Starting Matlab

Matlab for Windows Version 5 is used for all the examples in this handout. To start Matlab, select "programs" from the "start" button, go to Matlab menu set, select "Matlab" and click on it with left mouse button. To quit Matlab simply type »quit <return>. The symbol <return> is used throughout the text to denote pressing the return or enter key. Boldface characters denote both text to be typed and matrices. The context will indicate specific use.

Exercise 1: Start Matlab

Return to Top

2) The Command Line and Basic Commands

Once Matlab is started the Command Window appears as shown in figure. The command line with the prompt » appears inside this Window. Type your command and press Return. Matlab will interpret and execute your command.

The command line is the heart of using Matlab. All commands are entered at the command line, from adding two numbers, to executing programs (script files) and functions, or making graphs.

First consider how the command prompt can be used as if it was a hand calculator. To add two numbers, for example 5+7 enter the following at the command prompt:

»5+7 <return>

The answer pops out and is declared as the variable ans. If a variable name is not declared explicitly (such as c=5+7), Matlab always names the output ans. Hence, the latest results write over the previous contents of ans.

Next press the up arrow key and notice that the previous command comes up on the command line. Successive pressing of the up arrow displays earlier commands sequentially.

Now type »ans <return>. The number 12 should come up on the screen.

The basic commands for scalars are as follows:

    + addition
    - subtraction
    * multiplication
    ^ power
    / division

Entering calculations has the same precedence or order of evaluation as other languages. This means that commands are evaluated left to right with parenthesis with the highest precedence. For example lets say that you wanted to multiply 3 by 4 add 5 to that then divide by 6 + 4 and then divide all that by 2. This would be entered as follows:

»(3*4+5)/(6+4)/2 <return>

There are various other functions such as:

sin --> sine
cos --> cosine
tan --> tangent
exp --> exponential
log --> natural log
abs --> absolute value
There are many other functions available in Matlab. For example, to take the sin of pi/2 type

»sin(pi/2) <return>

where pi is already defined as 3.1416.

Next, instead of adding 5+7 lets two variables a+b where a=5, and b=7. First declare variables a and b, by typing

»a=5 <return>
»b=7 <return>
»a+b <return>

If you wish the variable c to equal a+b type »c=a+b <return>. More complicated expressions can be written; for example »c=(sin((a^2*b/tan(a)-b^2)/pi)*exp(a-b))/b. Matlab is case sensitive so the variables a and A, are considered as two different variables.

Type

»who <return>

A list of your variables in the workspace appears: a,b,c, and ans. The workspace represents which variables have been defined and are available for further calculations, manipulation or graphing. Now type »whos <return>. The command whos does the same as who, but in addition it will tell you how many elements are in each variable, for example the variables a,b,and c are scalars so there size is 1 by 1 with one element. If a was a column vector with 10 elements then its size would be 10 by 1.

The variables in the workspace can be saved to a file for future use. Matlab can save either all the variables in the workspace or selected variables. Before saving any variables to the disk type »pwd <return>. The command pwd (print working directory) will indicate the location (drive and directory) to which the variables will be saved. To save any variables to a floppy disk type »cd a: <return>. Next use the pwd command to check that the floppy drive a: has been properly set. Finally to save all the variables in the entire workspace type »save fl_test <return> where fl_test is the filename that the variables are stored and fl are your first and last initials. Note that between fl and test there is an underscore. The save command creates the file fl_test.mat in the current directory. If you like to save only certain variables just type the variables you wish to save after the filename, for example to save the variables a and b type »save fl_test2 a b <return>. If you desire to save data in ASCII then type »save fl_test3.dat a -ascii <return>. All data are stored contiguously in the ASCII version. You can not separate data for different variables in a file saved as ASCII.

To load files from disk use the load command. First empty the workspace of the old variables. Type »clear. The command clear empties the workspace of every variable. To clear a given variable or variables add the variable names after the word clear. Now to load the file name test type »load fl_test, (Make sure you are in the correct directory, use the pwd and cd commands if necessary). Next use the who command to check if the variables are in the workspace. To load an ASCII file type the complete filename after load including the extension ".dat". For example, »load fl_test3.dat.

You may make a copy of your session (all commands, messages from Matlab, data entered, and results) by activating the diary feature. Type diary filename to indicate the file where all session information will be stored as a text file, diary on to activate and diary off to terminate writing to the file. You may use this feature for your homework solutions by editing the raw information by using a word processor.

Return to Top

3) Script Files and Functions

Most of Matlab’s usefulness is in its wide variety of built-in functions. Some of these are in the computing engine of Matlab. For example sin and exp are such functions. Others are developed as files that can be called by a Matlab command with a statement of the form:

[outputs] = functionname(arguments)
where functionname is the function’s name and arguments are the inputs to the function. A function can have one or many arguments. The results (outputs) returned by a function are listed on the left side of the = sign, in square brackets ([…]). For example

» [V,D] = eig(X)

produces a diagonal matrix D of eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that X*V = V*D. There are many functions included in Matlab, please use the help system to learn the ones you may need. Matlab has an extensive help library to describe the use of functions. To use help use the mouse and the pull menu titled Help or simply type the command help followed by the function name. Typing help alone will give a list of subject areas. For example type »help <return>. A list of subject area should scroll down the screen. Look for the subject heading with the word general. Next type »help matlab\general <return>. Next a list of functions with a short description should follow. Next type »help load or another function name from the list. A detailed description then is scrolled down the screen.

A script file can be thought of as a program. It contains a series of commands or function calls. The difference between a function and a script file is that a function usually requires input arguments and will return only the output values. Just like a subroutine in Fortran or C, the variables used in a function are local. A script file requires any variables it needs defined within the workspace at the time of its use. Also a script file will return (to the workspace) every value it calculates. Script files are a simple way of writing programs and prototyping functions. They can be understood better through a simple example. The following steps will outline how to create a script file.

  1. From the Matlab command window point the mouse to the pull down menu File, then click on New, and finally click on M-file. A window with an editor should pop up. The screen should appear as follows (the actual window sizes vary and often they would overlap, this picture is arranged to show you both screens):

  2.  

     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     


     

  3. Move to the editor screen and type the following:

  4.  

     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     

    a=5
    b=7
    c=9;
    d=sqrt(a^2+b^2+c^2)
     

  5. Click on the pull down menu File of the Editor/Debugger window. Click on Save or Save-As. Under the Drives title choose the drive and folder (for example A:/) . Lastly, under File Name: write fl_first.m.

  6.  
  7. Go back to the Matlab window and type »fl_first <return>. DO NOT TYPE the extension ".m", otherwise you will get an error. The variables a,b and d with there corresponding values should scroll through the screen. If you get an error, probably the drive location is not set correctly. Use pwd to check the drive, and if not set on drive A use the cd command to go to drive A. Notice that variable c did not appear on the screen. This is because a semicolon follows the declaration of c in the script file fl_first.m. A semicolon will prevent any values assigned on that line from being displayed on the screen. If one of the variables was not declared in the script file, it would have to be declared in the workspace for the program to run. To investigate this, edit the file fl_first.m by removing the line b=7. Run fl_first.m again in the Matlab command window. You should get an error. Now type b=7 in the Matlab command window and press return. Run fl_first.m again. You should get the correct value of d on the screen.
Now develop a function to do the same task. It should have 3 inputs (x, y, z) and one output (q). When values or variables for x, y, and z are entered, the function should calculate q and rename it as the variable named in brackets. Go to the editor window and type:

function [q]= fl_firstfn(x, y, z)
q=sqrt(x^2+y^2+z^2)

and save the function as fl_firstfn.m. Go back to the Command window and enter a=5, b=7, and c=7 by using proper procedure. Now type:

» [d]=fl_firstfn(a, b, c)

The screen should show:

q =
    12.4499

d =
    12.4499

Exercise 2: Modify the function that you wrote so that q=12.4499 is not displayed when you execute the function.

To learn more about the function structure in Matlab, open the help window, type "function" in the upper left window (the current entry should be "MATLAB Help Topics") and press Enter.

Return to Top

4) Using Vectors and Matrices

Matlab is a powerful tool to manipulate arrays, vectors, and matrices. For example in C or FORTRAN77 the command to multiply two matrices would contain many lines of code and a "for" loop. In Matlab, the command to multiply two matrices A (for example a 2´ 3 matrix) and B (a 3´ 2 matrix) is simply A*B.

Define a 3´ 3 matrix A by typing »A=[1 2 3; 4 5 6; 7 8 9] <return>. The semicolon or a Return indicates the end the end of a matrix row. Now type »A(1,1) <return> . The first element of the first row and first column is displayed on the screen. Any single element of a matrix can be extracted in this manner. Next type »A(:,1) <return>. This extracts the first column. Any column can be extracted in this fashion. The same can be done for rows by replacing the column number by a colon (for example A(3,:) extracts the third row).

To extract a 2´ 2 matrix that are the elements of the first two rows and columns of A, type »A((1:2),(1:2)) <return>. The colon terms 1:2 mean 1 to 2 by step of size one, or 1,2. To understand this type »(1:10) <return>, and note what happens. Next type »(1:0.5:10) <return> and notice the effect of the step size 0.5 between 1 and 10.

Exercise 3: For matrices A and B defined below, type the following commands and observe the results:

»A=[5 6 7;10 4 6;2 8 10] <return>
»B=[1 3 0;5 13 5;8 2 1] <return>

»A+                                     »A-B
»A*B                                   »A^2
»B*A                                   »A*A
»5*A                                    »A.^2
»A’ Transpose of A        »A./5
»A^(-1)                                »A/A
»inv(A) Inverse of A       »A/B

Recall the "eig" function as an example of a matrix operation with an output argument:

»[V E] = eig(A) where V will be the eigenvectors of A and E its eigenvalues. Almost any function can be used on a matrix. For example, sin(A) would be a 2´ 2 matrix that has the value of sine of each element of A.

Enter the command »C=[A B]. C will be the 3´ 6 matrix that consists of A augmented with B. Why is this useful?

To learn more about matrix functions first type »help and press Return. A list of basic operators will scroll down the screen. For further information about matrix functions type »help matlab\elmat and »help matlab\matfun or go to the html of pdf help and electronic book files.

Exercise 4:

  1. Define a 3´ 3 identity matrix I by using the appropriate function. Hint: type »help matlab\elmat <return> and look for the function that corresponds to Identity matrix.

  2.  
  3. Define the following matrices:

  4.  
  5. A = [1 2 3;7 8 9;13 14 15] and B = [4 5 6;10 11 12;16 17 18]

  6.  
  7. Create the row vector C from A and B that is equal to C=[1 2 3 4 .... 18]. Hint: Start with C=A(?,?) then let C=[C B(?,?)] then C=[C A(?,?)],. Follow this pattern until you obtain the desired C. You fill in the correct values for "?".

  8.  
  9. What would be a really easy way to create the vector C without even using matrices A and B.

  10.  
  11. Create a 400´ 1 column vector of ones. Hint: Do not type [1;1;1;1;1;1;1...........]

  12.  
  13. Given matrix A from Part 2 find the mean values of the columns and rows using the mean function. Hint: Use help to learn how to use mean, and remember the transpose of a matrix.

  14.  
  15. Using A and B from Part 2 and D=[A B] determine which ones of the following will work:
a) A*D         g) A*B*D

b) D*A         h) A/B

c) A*B           i) A/D

d) D’*A         j) D/A

e) A*D’         k) D(:,(1:3))*A

f) A*D*B       l) D((1:3),:)*A
 

Matrix Operations (Section condensed from Georgia Papavasiliou’s Matlab Tutorial for Reaction Engineering)

There are two kinds of operations with arrays: matrix operations and operations on individual array elements. The first set includes taking the transpose of an array, ``matrix" multiply, adding and subtracting arrays of the same shape, etc. Matrix division is used to solve sets of linear equations. To solve the set of equations:

5x1 + 3x2 - x3 = 5
2x1 - 7x2 - 3x3 = 0
x1 + 5x2 + 6x3 = -7

that may be represented as: AX=B where A is a 3x3 matrix composed of the coefficients of x1, x2 and x3, and B is the column vector [5 0 -7]’.

We need to define a 3 by 3 matrix with the coefficients of the unknowns in it and divide a column vector with the right hand side coefficients by this matrix. There are two divide symbols in MATLAB.

A\B gives the solution to: A*X = B
A/B gives the solution to: X*A = B

In our case we need to use the first form:

» a=[5 3 -1 ; 2 -7 -3 ; 1 5 6];
» a\[5 0 -7]' <-- Note the quote or prime to transpose the vector.

ans =

    0.1168
    0.8426
    -1.8883

We can confirm our solution by using matrix multiplication:.

» a*ans

ans =
    5.0000
    0.0000
    -7.0000

Operations on individual array elements

 Let:

» a=[1 2 3];
» b=[2 5 8];

Suppose we try to multiply them:

» a*b

??? Error using ==> *
Inner matrix dimensions must agree.

What happened? The multiplication symbol by itself when applied to vectors or matrices is interpreted to mean matrix multiply. We could have done this by transposing the second vector:

» a*b'

ans =
    36

This gave (1*2 + 2*5 + 3*8) = 36, the usual scalar product of two vectors. If we want the element by element product, we must put a period before the multiplication symbol:

» a.*b

ans =
    2 10 24

The same procedure must be followed in doing division and exponentiation. For example, if you want the square of each element in a vector you must use:

» a.^2

ans =
    1 4 9

Return to Top

5) Plots and Graphics

Matlab can make two dimensional and three dimensional plots. Only the two dimensional plots will be covered in this tutorial. Two dimensional graphics utilize vectors and mainly the plot command. Type the following:

»plot(2:11,’*’) <return>

A graph is displayed with the y-axis values that go from 2 to 11. The x-axis goes from 1 to 10 reflecting the number of each element. Now try:

»plot(1:0.5:10,’*’) <return>

This should clarify further the values for the x-axis. Now try »plot(1:10) <return> and notice the type of line plotted. The term ’*’ in the previous plot command designates the type of line or point that are plotted. There are various line and point types such as ,’+’, ‘--’, ’-.’, ‘:’, ‘.’ etc. Use »help plot <return> for further details on line type.

Matlab can also plot two different values versus one another. For example, define »x=(1:0.1:100); <return> and find the sin of each value of x by »s=sin(x); <return>. To plot sin(x) versus x, type »plot(x,s) <return>. Now try »plot(x,sin(x)) <return>. Do you observe a difference? Now type »figure <return> and type »plot(sin(x)) <return>, what is the difference between the two figures? Lets find out by typing »hold on <return>, then »plot(x,sin(x)) <return>. Notice how the first plot appears on the second. The "hold on" command allows multiple plots on the same graph.

Close Figure No. 2 by using the mouse and the File pulldown menu (on the Figure No. 2 window). Use the xlabel, ylabel and title commands to label your graph. First type »xlabel(‘x from 1 to 100 by step 0.1’) <return> to label the x-axis, then type »ylabel(‘sin(x)’) <return> to label the y-axis, and lastly type »title(‘First Graph Sin(x)’) for adding the title.

There are many other useful graphing commands such as axis, subplot, contour, mesh, bar etc. Use the help command and specify either matlab\graph2d or matlab\graph3d. Here is a list of a few commands that may be helpful:

Function                 Description of Task Performed

clf                             Clears the graphics window
contour                     Allows you to generate a contour plot of a 3-D object
hold                          Holds the current plot so you can add to it
subplot                     Allows you to put many in the same window
figure                        Allows you to create multiple graphics windows
close                         Allows you to close graphics windows
plotyy                       Creates 2d plot with different y axes on 2 sides
loglog                       Creates a log-log plot
semilogy/semilogx     Create semi-log plots with one coordinate on a log scale.
gtext                         Place text in a plot with the mouse
ginput                       Get data using the mouse from a graph
 

To save your graph to disk use the following command sequence: »pwd <return> --> Make sure your in drive a: if not type »cd a: <return>
»print fl_grph -dmfile <return>
Note that you must specify your file name fl_grph, type a space, and then type -dmfile.

The print command will save to files to your floppy disk, fl_grph.m and fl_grph.mat. These files are needed to reproduce the graph so do not delete them. Next close the figure window, the window with the graph and type » fl_grph <return>. The graph should come up to the screen as before.

It will be important that you know how to save variables, graphs and text files, and generate dairy files so that you can turn in your homework on a floppy disk or as e-mail attachment.

Exercise 5: Define a vector x that goes from 0 to 100 by steps of size 0.5. Plot sin(x) and 1+cos(x) on the same graph. Represent the sine curve with a ‘dashed’ line and the cosine curve by a ‘dashdot’. Also label the graph and save it as ‘fl_hwe1’.

Exercise 6: Graph a ramp function. A ramp function consists of a constant slope that ramps or slopes upward or downward over time. After a specific ramping period, the function may reach a final value and remain at that value for a period of time. A physical example would be the temperature of water that is being heated. It will increase until it starts to boil then remain constant as long as water and steam coexist. If the system is heated further all the water may evaporate to yield eventually superheated steam (temperature will start ramping again after all the water has evaporated).

Plot a ramp that starts at time 10 at a value of 5 and ends at time 20 with a constant value of 25. Hint: use the zeros and ones functions. Create a vector which has X values and resembles [5 5 5 . . . . . .7 9 . . . . 25 25 25 ].

Return to Top

6) Programming in Matlab

When writing a script file, a function, or essentially a program, Matlab has the familiar programming tools such as ‘for’ loops, and ‘if’ and ‘while’ statements. The basic structure of each is as follows:

The ‘if’ statements

if ‘variable’ ‘conditional statement’ constant
          (statements)
      elseif ‘variable’ ‘conditional statement’ constant
           (statements)
       else
           (statements)
end

the if statement only requires the first ‘if’ and the end statement. The ‘elseif’ and ‘else’ can be added if additional checks are required.

The ‘for’ loop

for ‘counter’
    (statements)
end

The ‘while’ loop

while ‘conditional statement’
    (statements)
end

A simple example:

for I=1:10
    if I <= 5
        c(I)=I+1
    elseif I < 7
        c=[c I]
    else
        c=c.^2
    end
end

Another example of a ‘for’ and ‘while loop combination is as follows:

TEST=zeros(3,3)
k=0
for I=1:10
    while k = = 0
        TEST=TEST+eye(3);
        if trace(TEST) < 20
            k=0;
        else
            k=1;
        end
    end
end

A semicolon should not be placed after for, while, if, and end statements. The semicolon should only be placed after manipulations or assignments of variables. Also notice above how the variables TEST and k are predefined before the for loop and while statements.

Exercise 7: Write a program that creates the ramp function in Exercise 6 by using a for loop and if statements.

Return to Top

7) Solving Ordinary Differential Equations

Matlab has several functions for numerical solution of ordinary differential equations:

ODE23- Solve non-stiff differential equations, low order method. 2nd/3rd order Runge-Kutta method (medium accuracy).
ODE45- Solve non-stiff differential equations, medium order method. 4th/5th order Runge-Kutta method (high accuracy).
ODE23S- Solve stiff differential equations, low order method.
ODE15S- Solve stiff differential equations, variable order method.

Examples:

Example 1. Consider the following system of first-order ODE’s:

The solution to Matlab requires two programs to be written. The main program contains the initial conditions, time of integration, and a call statement to the second program that contains the differential equations. The first program is entered. After starting Matlab by double clicking on the Matlab icon select FILE/NEW/M-FILE in the Matlab command window to create the file containing the differential equations. Enter the following lines: function DXDT = fl_test(t,x)
dx1dt = x(1)*(1-x(2)^2)-x(2);
dx2dt = x(1);
DXDT = [dx1dt dx2dt]’;
Save this file by clicking on FILE/SAVE-AS under the name fl_test.m. If you have already made a file with that name, use fl_test1.m or some similar file name. Next, open a new m-file and enter the following code: x0 = [0 0.25];
TSPAN = [0 20];
[t,x] = ode23(‘fl_test’,TSPAN,x0);
plot(t,x(:,1),’ *’,t,x(:,2),’o’)
title(‘x vs. Time’)
ylabel(‘x’)
xlabel(‘time’)
legend(‘x1’,’x2’,0)
Save this file as fl_testsim.m. Run the programs by entering fl_testsim (without the .m extension) in the Matlab command window. The program will generate the following figure.
   Example 2. Consider the irreversible, exothermic reaction  in a CSTR as outlined by Uppal, Ray, and Poore:

where x1 is the conversion of the reactor and x2 the dimensionless reactor temperature. Now to solve these equations using ode23 we must first represent these ODE in an m-file. To define the CSTR open a new m-file, call it cstr.m and type the following in your m-file:

%Exothermic, irreversible CSTR A --> B ODE's
%Define your function
function DXDT = cstr(t,x)
%Set constants, note these could also be algebraic equations such that
%they are functions of your states x1 and x2
B=22.0;
beta=3.0;
gamma=inf;
x2c=0;
Da=0.085;
%declare ODE's
dx1dt = -x(1) + Da*(1-x(1))*exp(x(2) / (1 + x(2) / gamma));
dx2dt= - x(2) + B*Da*(1-x(1))*exp(x(2) / (1+x(2) / gamma)) - beta*(x(2) - x2c);
%Make vector of derivatives
DXDT=[dx1dt dx2dt]';

Note DXDT is a vector containing the derivatives of the states while the vector x are the values of the state variables themselves. Next open a new m-file and call this one cstrsim.m. In this file we will declare the integration region T0 and Tfinal, and the initial conditions X0. Type the following (program developed by Jeff DeCicco):

%File for initiating simulation of CSTR and plotting time responses and
%Phase Plane Plots
%Declare initial conditions and time interval for integration
Xo=[0.9996 10.8]
To=0;
Tfinal=20;
TSPAN= [To Tfinal];
%Solve ODE's using ode23
[T X]=ode23('cstr',TSPAN,Xo);
%Plot yout output
figure
subplot(2,1,1)
plot(T,X(:,1))
xlabel('Time')
ylabel('Conversion')
title('Dimensionless Conversion x1 versus Time t')
subplot(2,1,2)
plot(T,X(:,2))
xlabel('Time')
ylabel('Dimensionless Temperature')
title('Dimensionless Temperature x2 versus Time')
figure
plot(X(:,1),X(:,2),'g')
hold on
%Solve same ODE for different initial conditions
Xo=[0.93 5.1]
[T X]=ode23('cstr',TSPAN,Xo);
plot(X(:,1),X(:,2),'y--')
axis([0.7 1.1 3 9])
title('Phase Plane Plot Fig. 10 of Uppal, Ray and Poore')
xlabel('X1')
ylabel('X2')
figure
subplot(2,1,1)
plot(T,X(:,1),'--')
xlabel('Time')
ylabel('Conversion')
title('Dimensionless Conversion x1 versus Time t')
subplot(2,1,2)
plot(T,X(:,2),'--')
xlabel('Time')
ylabel('Dimensionless Temperature')
title('Dimensionless Temperature x2 versus Time')

To simulate the program type »cstrsim <return>. It may take a few minutes for the simulation to complete its calculations. Time is in vector T, and conversion and temperature are in the n´ 2 matrix X. Thus you can plot the state variables by typing plot(T,X(:,1)) for conversion or plot(T,X(:,2)) for temperature as shown in the program. The plots generated are:


 

The function ode45 will lead to more accurate results but may take a significantly longer time to finish its calculations than ode23.

Return to Top
 

8) Process Control

Matlab has several toolboxes with functions that facilitate process model identification, system stability assessment, control system design, and process simulation. The toolboxes deal with a wide variety of control areas from conventional control to optimal and model predictive control. The most frequently used toolbox in the undergraduate process control course is the Control System toolbox. In graduate courses, the System Identification and Model Predictive Control toolboxes are also used.  To see a list of the functions available in the the Control System toolbox type »help matlab\control <return> or use the other help options. These functions primarily work with model representations that are either transfer functions or state space. The following will give a brief outline of some common tasks in defining a model and using some common functions.

First consider a simple SISO second order transfer function in the s domain:

Basically the transfer function is the ratio of two polynomials which can conveniently be represented by their coefficients in decreasing order in MATLAB. Label the numerator polynomial as num and the denominator polynomial as den by typing:

»num=[1 1] <return>
»den=[4 2 1] <return>

Now num and den represent our model.  The convention for declaring a polynomial is that the first element of the array is the coefficient of the highest order in s and the last element of the array is the constant term. Thus for num that has two elements the highest order is 1 because the last number will always represent a constant. For den the last number still represents the constant so the highest order is two. Most Matlab functions can operate on a transfer function that can be defined by using the function tf:

»g=tf(num,den) <return>

One can also define the transfer function g by writing directly:

»g=tf([1 1], [4 2 1]) <return>

As an exercise type:

»num2=[1] <return>
»den2=[4 2 0] <return>

What would be the structure of this transfer function? Remember the last number is always a constant and then count from right to left in increasing order in s. Check by using tf that the transfer function is:

A task that would be of interest is to check the stability of the system by finding the poles of the transfer function, and the system response to certain inputs. To check the poles use the roots function. Type »help roots <return> and read the description. Next type »roots(den) <return> followed by »roots(den1) <return>. Now you should be able to describe the behavior of each transfer function by looking at the poles.  The same task can be achieved by the command »pole(g).

Now consider the state space representation:

where x is a  vector, A matrix, B matrix, C an  matrix, and D an  matrix. Thus k denotes the number of state variables, m the number of outputs, and n the number of inputs. To define a state space model the matrices A, B, C, and D must be specified. Every matrix has to be defined even if a matrix is a zero matrix. How do you check the stability of a state space model? Try »help eig <return>.

Now lets investigate the power of the control toolbox. Type »help toolbox\control <return>. Take a short look at the list of functions. One of the most important tasks is constructing the open-loop or closed-loop transfer functions from the models of the process and control loop equipment. There are various functions in the control toolbox that allows one to combine two transfer functions, or close a loop, or add new dynamics. No one function can construct a closed-loop transfer function from many different individual transfer functions. To do this one must use a combination of several functions such as, series, cloop, etc. All of these functions are listed under the subtitle Model Building after you type »help toolbox\control <return>. They are fairly straightforward to use. However, to make life easier the Control System toolbox has m-files append and connect that will build a state space representation of the closed-loop system.  A detailed example for developing the closed-loop system is given in following section.

Given one type of model, one may be interested in converting it to another form. Control System toolbox has functions to convert most models to each other. For example, the function ss2tf converts a state space model to a transfer function model. Other model conversion functions are listed under Model Conversions subtitle in the help file.

Control System toolbox has functions to compute the transient response of an open or closed loop process to various input changes. The functions step and impulse compute and plot the responses to step and impulse changes. A more general function, lsim, computes the response to any input defined by the user.

Example:  Compute of the poles of G(s)=(3s+1)/(s^2+2s+1) and the response of the system described by G(s) to a step change of magnitude 2. Plot the response.

% Computation of the poles of G(s)=(3s+1)/(s^2+2s+1)
% Enter the denominator
den=[1 2 1];
% Find roots of denominator.  They are the poles of G(s)
poles=roots(den)

% You may also enter the transfer function and compute the poles directly
g1=tf([3 1], [1 2 1])
poles2=pole(g1)

% Computation of the response to a step change of magnitude 2
% Assume that the values of the response will be computed for a total time of
% 20, in increments of 0.1.  Generate the time vector
t=0:0.1:20;

% Method 1: Use the function "step".  If there is no left side argument such as
% y=step(transfer function, time), then Matlab generates the plot automatically.
% If you wish to store the output values and/or give specific headings in your
% graph, include the left side argument then make the plot by including appropriate
% commands. Then the commands would be:
%  y=step(g1,t)
%  xlabel('time')
%  ylabel('Output')
%  title('Plot of output as a function of time to a step input of magnitude 2')
%  plot(t,y)
step(g1,t);
input('After inspecting the figure, close it and press ENTER to continue');
pause

% Clear the current figure
clf
% Method 2:  Use the general simulation function "lsim".  You must specify
% the input sequence as a vector of the same length as the time steps
% Determine the number of elements in t by using the command size(t).
% Label the input sequence as u:
u=2*ones(size(t));
y=lsim(g1,u,t);
y=step(g1,t);
plot(t,y)
xlabel('time')
ylabel('Output')
title('Plot of output as a function of time to a step input of magnitude 2')
% If you use lsim without the left side argument, Matlab generates the plot automatically

It is important to know if a system is stable or unstable.  This can be done by using various numerical and graphical methods.  The numerical methods are based on the computation of the poles or eigenvalues of the system.  If the real part of a pole or eigenvalue is positive, the system is unstable.  Three graphical procedures that are used for closed-loop SISO systems are root locus, Bode, and Nyquist plots. The corresponding functions are rlocus, bode,and nyquist.  There is a function that can report the controller gain that corresponds to a specific root indicated on a root locus plot, rlocfind.  It can be used only after a root locus has been plotted by an rlocus command.  Note that the open-loop transfer function or state space model is written in the arguments of these functions, but the resulting plots are for the closed-loop system.

Return to Top

9) Closed-loop Model Development Example: Mixing Process Model

The example illustrates the use of the functions append and connect in the Control System toolbox to construct the open loop and closed-loop models for a mixing process. A schematic of the process is given in Figure 1 (This example is based on Problem 10.15 in Seborg, Edgar and Mellichamp, 1989) and its block diagram is shown in Figure 2.
 


Figure 1.  Mixing Process Problem
 

Figure 2.  Block diagram of mixing process and control system

The goal is to build the state-space model of the overall (closed-loop) system. The models of various parts of the process may be in different forms as illustrated in the help files for connect. The transfer functions representing the mixing process are:

To define our system, open a new m-file and save it as fl_mix.m where fl are your first and last initials. Declare all the variables and following commands in your m-file.  The transfer function for the time delay can not be directly represented in MATLAB. The transfer line represents a time delay in the measurement of concentration. Recall MATLAB accepts only ratios of polynomials or state space representations, so we will have to use a Pade approximation of the time delay. Type »help pade <return> and read the description. To define the transfer line use [n8 d8]=pade(0.52,2). This represents a second order Pade approximation (the numerator and denominator are second order) with time delay of 0.52 minutes. To keep track of system components you may assign a number to each block and use it in your program statements.

Now define the system.. Type the following in your m-file.

n1=0.08; % the setpoint transfer function Km
d1=1;
sys1=tf(n1,d1);
n2=[2.016*0.87*0.154  2.016*0.87  2.016];
d2=[0.87   0];
sys2=tf(n2,d2);
n3=0.75;
d3=1;
sys3=tf(n3,d3);
n4=0.082;
d4=[0.2 1];
sys4=tf(n4,d4);
n5=93.3;
d5=[0.67 1];
sys5=tf(n5,d5);
n6=0.93;
d6=[0.67 1];
sys6=tf(n6,d6);
n7=1;
d7=1;
sys7=tf(n7,d7);
[n8 d8]=pade(0.52,2);
sys8=tf(n8,d8);
n9=0.08; % Transfer function for Gm
d9=1;
sys9=tf(n9,d9);

Notice that there has been a block added n7, d7 that is simply an all-pass transfer function i.e. 1/1 on the schematic for our process. The reason for adding this block will be clear when function connect is explained.

To have a physically realizable transfer function the order of its numerator has to be less than or equal to the order of its denominator. We can not create the transfer function for the controller as we wrote in n2 and d2, since it is not physically realizable. Basically our PID controller is physically unrealizable and we have to approximate it with

where is some small number between 0.05 and 0.02. So we must redefine our PID controller as follows:

n2=[2.016*0.87*0.154 2.016*(0.87+0.154) 2.016]
d2=[0.05*0.87*0.154 0.87 0]

For practice you may enter instead the transfer functions of the approximate PID controller and let Matlab do the controller model generation for you automatically.

Combine different system parts by using function append.

sys=append(sys1,sys2,sys3,sys4,sys5,sys6,sys7,sys8,sys9);

Basically append constructs a state space representation given in a,b,c,d. Note that a,b,c,d does not represent our open or closed loop system. To construct our closed loop system the function connect must be utilized.  Type »help connect <return> to learn more about this commandThe command append has produced the matrices a,b,c,d, now we have to specify the connections between various blocks by defining and array Q, and the inputs and outputs to the system. The variable Q simply states how our system blocks are interconnected. From the above description you should be able to understand that Q will defined as follows.

Q=[2 1 -9
3 2 0
4 3 0
5 4 0
7 5 6
8 7 0
9 8 0]

The first line 2 1 -9 indicates that block 2 (the PID controller) has an input which is the sum of the output from block 1 minus the output from block 9. The line 3 2 0 indicates that block 3 has an input from block 2. The other rows are also developed using the block diagram connections until all the loop elements are accounted for.

Next we need to specify which blocks are the inputs and outputs of our system. Type

inputs=[1 6]
outputs=7

Now it should be apparent why block 7 was added. We need a reference block so that we can extract the output. Blocks 5 or 6 could not be our output because our output is the sum of the outputs from these blocks and that is exactly what block 7 is defined to be. Lastly our inputs simply come from blocks 1 (setpoint) and 6 (disturbance).

We have everything needed to run connect, so type in your m-file after inputs and outputs:

sysc = connect(sys,Q,inputs,outputs);
sysr = minreal(sysc);

sysc contains the closed-loop system model. The last function, minreal was added since the function connect does not necessarily give the minimum realization state-space model.

We have finished the m-file for developing the closed-loop model of the mixing process, save your m-file and then run it. Then experiment with functions such as bode, step, impulse, etc. Remember that you have a state space model with two inputs (one for the reference signal and one for the disturbance).  For example, the command

step(sysr)

gives the following figure that shows the response of the output to a unit step change in the setpoint (reference) and a unit step change in the disturbance.


 

Return to Top

References:

Seborg, D., T. Edgar, and D. Mellichamp, Process Dynamics and Control, Wiley, 1989.

Uppal, A., Ray W. H., and Poore A. B., ‘On the Dynamic Behavior of Continuous Stirred Tank Reactors’, Chem. Engng Sci. 29, 967-985, 1974.