Master Scilab!by 2. April 2008 Contents:
1 What is Scilab? 1 What is Scilab?Quoted from the homepage of Scilab at http://scilab.org: Scilab is a free scientific software package for numerical computations providing a powerful open computing environment for engineering and scientific applications. Scilab is an open source software. Since 1994 it has been distributed freely along with the source code via the Internet. It is currently used in educational and industrial environments around the world. Scilab includes hundreds of mathematical functions with the possibility to add interactively programs from various languages (C, C++, Fortran…). It has sophisticated data structures (including lists, polynomials, rational functions, linear systems...), an interpreter and a high level programming language. Scilab is quite similar to Matlab, and the range of functions are comparable. The largest benefit of Scilab is of course that it is free :-). Also, Scilab is easy and fast to install (and you do not have to restart your PC before starting to use it). Scilab is also similar to Octave, which is also free! Octave is more similar to Matlab than to Scilab. One problem with Octave has been that data plotting is more cumbersome in Octave than in Scilab. (You can have both Scilab and Octave installed :-) One nice thing about Scilab is that you get Scicos automatically installed when you install Scilab. Scicos is a block-diagram based simulation tool similar to Simulink and LabVIEW Simulation Module. 2 About this documentThis tutorial guides you through the steps towards mastering Scilab. I have written this document because I did not find a proper tutorial on the Scilab homepage. I assume that you do all the activities in the blue boxes, as here:
Please send comments or suggestions to improve this tutorial via e-mail to finn.haugen@hit.no. 3 Downloading and installing ScilabThe installation file, which is an *.exe-file, is available for download at http://scilab.org. Once you have downloaded this exe-file, open (run) it, and then follow the instructions on the screen. (It should not be necessary to restart your PC before you start Scilab after the installation.) Note that by installing Scilab, you also get Scicos installed. 4 The Scilab environmentTo start Scilab:
Starting Scilab opens the Scilab command window, see the figure below.
The Scilab command window Scilab commands are executed at the command line by entering the command, and then clicking the Enter button on the keyboard.
The result is shown in the command window (see the figure above). 5 Scilab Help
The Help window is shown below.
Scilab Help window As you see from the Help window, the commands and functions are organized in a number of categories.
The functions are as shown in the figure above. To get detailed help text about a specific function, click that function.
The detailed help text for the abs function is shown in the figure below.
The detailed help text for the abs function You can also search for a function by first clicking the Search button in the Help window (the magnifying glass button).
The result of the search is a list of relevant functions, see the figure below.
The result of the search for sine 5 Basic Scilab operationsTypically you use variables in your calculations. To create the variable a and assigning to it the result of 1+1:
Hereafter, (Enter) will not be shown, but it is assumed that you click the Enter button. The response is shown in the command window (but shown here). Now, try (remember to type the semicolon):
The response is not shown in the command window. The command was actually executed, but because of the semicolon the response was not shown. To verify that the variable b actually exists:
As you will see, it exists. If you omit the variable name, the result is assigned to the inbuilt variable ans:
The response shows that the variable ans has got the value 4. You can execute one or more commands - separated by semicolon or comma - on one line:
Scilab is case-sensitive. For example, d and D are two different variables:
As you see from the response (not shown here), d exists, while D does not exist (since we have not created D). Scilab variables exists in the workspace. There are two ways to see the contents of a workspace:
The response should be similar to what is shown in the figure below. (The user-defined variables are shown among many other variables.)
The response of the command who
This opens the Browser Variables window, see the figure below.
Browser Variables window The Browser Variables window contains at the bottom a number of utility buttons (not described in detail here). Note that if you exit from Scilab, the variables you created in the workspace are deleted. You can save variables in a file using the save function. However, if you really need to save variables that are a result of some Scilab expressions, then you should consider writing these expressions in a Scilab script instead. More about scripts soon. There are various ways to enter numbers (the pi is an inbuilt constant). Here are some illustrative examples (I assume that you see the principles from these examples):
The response is shown in the figure below.
Various ways to enter numbers You can determine how numbers are displayed in the command window with the format function, but the internal representation of the number in Scilab is independent if the display format. We will not look at details. If you need to change the display format, consult Scilab Help. Scilab functions are vectorized, i.e. functions can be called with vectorial arguments. (A vector is simply a one-dimensional matrix. We will return to vector- and matrix operations in a later section.) In the following example, first a vector of name t is created, then this vector is used as an argument in the sine function (the sine function assumes the argument is an angle in radians).
The response is shown in the figure below.
The result of the vectorized function call sin(0.1*t) where t is a vector 6 ScriptsA Scilab script is a text file of name *.sce containing Scilab commands. You can edit the script using the inbuilt Scipad editor. (Scripts can also have names *.sci. The default name when saving a fle in Scipad is *.sce.) You should use scripts even for small tasks because in this way you have all your "projects" saved in files which is good for documentation and also very convenient when you want to run all your commands after some changes. We will now create a simple script, and then run it. Running a script is the same as executing all the commands (from top to bottom in the script) at the command line one by one.
The Scipad editor is shown in the figure below. Note that double slashes (//) are used to start comments in the script.
Scilab script of name script1.sce opened in the Scipad editor Note that you can open several scripts in the same Scipad window with the File / New menu.
There are two ways to run the script1.sce script:
Let us try the Execute menu first:
The result is shown in the command window. Now let us try the exec command:
Scilab responds with an error! See below.
Scilab responds with an error when after executing the command exec script1.sce The error comes because script1.sce is not in the Current directory of Scilab. The Current directory is the directory where Scilab looks for your script when you try to execute it with the exec command. What is then the present Current directory?
The response (in the command window) may be different on different PCs. On my PC the response is
We saved script1.sce in the C:\temp directory which is different from the present Current Directory. Let us change the Current Directory to C:\temp, and then execute the script:
Now the script should run without errors. 7 Matrix operationsIn Scilab numbers are in general stored in matrices, which can be regarded as a table. Matrices with only one column or only one row are denoted vectors. (Matrices and vectors are very similar to arrays as used in programming languages as C, Visual Basic etc.) Creating matrices. Retrieving data from matricesLet us create a 2x2 matrix of name A. In general, comma (or space) is used to separate numers on a row (or line), and semicolon is used to separate rows.
The figure below shows A as displayed in the command window.
Creating and displaying the matrix A Create a row vector with first element 0, last element 1 and increment (step size) 0.2:
Create a column vector, note the apostrophe to transpose:
Create matrix B from given column vectors x and y:
Get the size (dimensions) of matrix B:
Get the first element of vector r. Note that 1 (not 0) is the first index!
Now try
You get an error because index number 0 can not be used. Get the second column in matrix B, and assign the values to vector v:
Some special matricesCreate an identity matrix of dimension 2x2:
Create a matrix of dimension 3x2 consisting of ones:
Create a matrix of dimension 3x2 consisting of zeros:
Matrix calculationsMatrices can be used in matrix calculations. Add two matrices (assuming A and C have been created as explained above):
Multiply two matrices:
Take the inverse of a matrix:
Elementwise calculationsElementwise calculations are made with the dot operator. As an example, to perform elemenwise multiplication of two vectors:
The result is shown in the figure below. Element z(1) is v1(1)*v2(1), and element z(2) is v1(2)*v2(2).
Elementwise multiplication of two vectors 8 PlottingThe basic command for plotting data in a figure is plot. In the following are a number of examples demonstrating various plotting options. (It is of course up to you if you type the code from scratch on the command line or in Scipad, or if you just copy the code. If you type the commands, you may omit the comments which is the text following the // sign to save time.) First, let us generate some data that will be used in the following examples.
Very basic plotting:
Below is shown the Scilab figure with the plot. Along the x-axis are the indexes of the y vector. The indexes are integers from 1 to 101.
Before we continue with more plotting commands, let us take a look at some buttons and menus in the Graphics window.
This opens the Clicking the GED button opens the Graphics Editor, see the figure below.
The Graphics Editor With the graphics editor you can change line colours, line style, add labels to the axis, add grid, etc. The various options will not be described here because it is rather easy to investigate the possibilities by yourself. Many of the options in the Graphics Editor can alternatively be set with options to the plot command. This will be shown in subsequent examples. You can produce various graphics files from the plot:
This opens the Export dialog window shown below.
The Export dialog in the figure window If you want to create a graphis file to put into a document processor, as MS Word or Scientific Workplace, you should select Enhanced Meta File (EMF), whch is a vectorized graphics format which means that the picture can be enlarged and still look sharp. However, EMF files can not be used in native web documents, e.g. in HTML-files to be displayed in a web browser. In this case you should select the GIF format (this format does not give you vectorized graphics). We continue with looking at more options to the plot command. Assume that we will plot y against t in Figure 1 which is the same figure as we used above. This is done with the command plot(t,y) where it is of course assumed that vectors t and y have the same length (same number of elements). If you just use the plot command, the new plot adds to the previous plot, showing two (or more curves). Typically, this is not what you want. To clear the previous plot, we use the clf (clear figure) command before we use the plot command.
The result is shown in the figure below. Observe that the x-axis now contains the t values.
Suppose you want to show the plot in a new Figure 2 in stead of the previously opened Figure 1:
To clear Figure 2:
To plot with with grid and plot title and labels:
The resulting plot is shown in the figure below. (If you think the code ax1=gca();ax1.grid=[0,0]; is too cumbersome for generating a grid, you can in stead use the Graphics Editor (by clicking the GED button in the figure window).)
To plot two curves in one plot (note the ordering of the vectors in the plot command):
The result is shown in the figure below.
To plot two curves with options for the line style and colours:
(Many more options can be found from Help linespec.) The result is shown in the figure below.
To plot several subplots in one figure:
The result is shown in the figure below.
Subplots 9 Functions for dynamics and controlIn the following are examples of using Scilab functions for solving problems in the field of control engineering. The functions used are among the function in the General System and Control function category (see Scilab Help). 9.1 Simulation of continuous-time transfer functionsThe csim function is used to simulate continuous-time transfer functions. Here is one example:
The two plots are shown in the figures below.
Step response
Sinusoid response 9.2 Frequency response of continuous-time transfer functionsThe bode function plots frequency response in a Bode plot. Here is one example:
The Bode plot is shown in the figure below.
Bode plot The bode function does not accept models with time-delay. However, there is a workaround to this problem that can be downloaded from here. 9.3 Simulation of discrete-time transfer functionsThe flts function can be used to simulate a discrete-time transfer function. Here is one example:
Comment to the transfer function in this example: It is a*z/(z-(1-a)) = y(z)/u(z) where y is output and u is input. This transfer function is "by chance" the transfer function corresponding to the difference equation y(k) = (1-a)*y(k) + a*u(k) which is the well-known EWMA lowpass filter (Exponentially Weighted Moving Average) frequently used in applications where a simple filter is needed, as in control systems where a lowpass filter is used to attenuate measurement noise from the sensor. The simulated response is shown in the figure below.
Step response 9.4 Frequency response of discrete-time transfer functionsThe bode function plots frequency response in a Bode plot. Here is one example:
The Bode plots are shown in the figure below.
Bode plots of discrete-time transfer function Comments about frequency information: The frequency axis shows normalized frequencies. The red vertical line in the Bode plots represents the normalized Nyquist frequency, which is 0.5, corresponding to the normalized sampling frequency being 1 samples/sec. If the actual sampling frequency is f_{s} [samples/sec], the actual (non-normalized) frequency is found by multiplying the normalized frequency by f_{s}. 9.5 Simulation of continuous-time state-space modelsThe state-space model is assumed to be dx/dt = A*x + B*u y = C*x + D*u with x0 as initial state. In the following example, first the state-space model is defined, then the responses due to initial state and to external input signal u are simulated using the csim function, and finally the responses in output y and in state-variables x are plotted.
The figures below show the response in the output y and in the state-variables x, respectively.
The response in the output y
The response in the state-variables x 9.6 Discretizing continuous-time systemsThe function dscr is used to discretize continuous-time models. dscr only handles state-space models. So, if you are to discretize transfer functions, you have to convert the original continuous-time transfer function into a state-space model (using the tf2ss function) before you use dscr. And thereafter, you convert the discrete-time state-space model to the (final) transfer function (using the ss2tf function). Here is one example:
The original continuous-time transfer function TFcont and the resulting discrete-time transfer function TFdisc are shown in the figure below.
The original TFcont and the resulting TFdisc Although it is not written in the help text of dscr you can assume that the discretization is based on the exact discretization with zero-order hold (piecewise constant input during the time steps) on the input of the system. (The result of using the c2d function in Matlab with zero-order hold one the above example gives the same result as above.) 9.7 Deriving transfer functions from state-space modelsThe function ss2tf derives transfer function(s) from a state-space model:
The resulting transfer function is shown in the figure below.
The resulting transfer function TF 9.8 Combining models: Series, parallel, and feedbackHere are examples of how to calculate the combined transfer function for series, parallel, and feedback combination of models (here transfer functions):
The results are shown in the figure below (you should verify the results by hand calculations).
Results of series, parallel, and feedback combination of the transfer functions sys1 and sys2 9.9 Frequency response analysis and simulation of feedback (control) systemsThe gain margin and phase margin and the corresponding frequencies can be calculated with g_margin and p_margin functions. Here is an example:
The calculated stability margins and crossover frequencies are shown in the figure below. The gain margin and the phase margin are larger than zero, and hence the control system is asymptotically stable.
Calculated stability margins and crossover frequencies The Bode plot of the open-loop system is shown in the figure below.
Bode plot of open-loop system The process output and its reference is shown in the figure below.
Simulation of closed-loop system The eigenvalues of the closed-loop system are calculated with the spec function, cf. the script above. The result is shown in the figure below. All eigenvalues have negative real part, and hence it is confirmed that the control system is asymptotically stable.
Eigenvalues (poles) of the closed-loop system 9.10 LQ (linear quadratic) optimal controllerGiven the following discrete-time state-space model: x(k+1) = A*x(k) + B*u(k) The lqr function (linear quadratic regulator) can be used to calculate the steady-state controller gain, K, of the controller u(k) = K*x(k) that minimizes the linear quadratic cost function J = sum_{0}^{inf}[x'(k)*Q*x(k) + u'(k)*R*u(k)) There is one problem with the lqr function, namely that the user can not define the weight matrices Q and R directly as arguments to the function, but in stead some intermediate calculations have to be done. A more user-friendly way to solve the control design task is to use the Riccatti solver ricc in Scilab [reference] (an Riccatti equation is solved as a part of the design algorithm). Here is one example:
The controller gain and the eigenvalues are shown in the figure below.
Controller gain and eigenvalues of LQ optimal control system (The Matlab function dlqr gave the same result.) 9.11 Kalman Filter gainsSteady-state Kalman Filter gains for continuous-time systems and for discrete-time systems can be computed with the lqe function. For discrete-time systems the Kalman Filter is assumed to have the one-step form (combined prediction and measurement update), and not on the probably more popular two-step form (separate measurement update and prediction update). |