Introduction to Matlab/Simulink

Department of Machine Design

HT2022
Preface

This text began as a handful of pages, highlighting a few useful commands for system
modelling and control design in Matlab. Since then, it has grown several times its original
size, with several additional chapters on advanced topics, and the language gradually
taking on a tone suitable for a comprehensive tutorial.
Despite this, the purpose of this document has remained the same – to serve as a handbook
to getting started with Matlab, understanding core concepts, producing high-quality fig-
ures, using the numeric LTI models included in the Control System Toolbox, and graphical
programming in Simulink. The idea is to serve as a tutorial on using the above mentioned
tools within the scope of simulation and control of dynamic systems. Small stand-alone
exercises are integrated into the reading material to assist the learning process.
This text has been tailored to support the courses given by the Department of Machine
Design at KTH Royal Institute of Technology, but will likely be useful in a larger scope, and
can be distributed freely without the author’s consent.
Prerequisites: This is not meant to be an introduction to programming, it is assumed
that the reader has knowledge in basic programming and computer science, numerical
methods, and algebra. For the chapters on dynamic systems modelling, basic knowledge in
automatic control (sometimes designated “signals and systems”) is strongly recommended.
Having worked in Matlab previously is not a requirement, but will be beneficial, in which
case the first chapters may be repetition.

ii
Acknowledgements

Some sections of this text are based on the tutorials Einführung in MATLAB and Einführung in
Simulink, originally written for the course Regelungstechnisches Labor, held by the Institute
of Automatic Control (IRT) at RWTH Aachen University in Aachen, Germany during the
summer of 2017. The author of this document would like to thank Prof. Dirk Abel and his
staff at the IRT for providing such excellent tutorials in the German language.
The opening sections on Simulink, including the material on running a demo model are
adapted from an older text from the Department of Machine Design.
Finally, the author would like to express heartfelt gratitude to all the people who helped
test the material, and gave their feedback and support.

iii
Contents

1 Introduction 1

2 Fundamentals 2
2.1 Introducing the desktop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2 The command line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.1 Suppressing output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.2 Input history . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.3 Auto-completion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2.4 Documentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2.5 Aborting commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2.6 Unquoted strings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2.7 Housekeeping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 Scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3.1 Creating new scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3.2 Sectioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3.3 Debugging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3 Basic functions and datatypes 7

3.1 System variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.2 Elementary mathematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.3 Vectors and matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.3.1 Declaring matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.3.2 Diagonal matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3.3 Sparse matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3.4 Accessing elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3.5 Matrix operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.3.6 Element-wise operations . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3.7 Matrix functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.4 Logical values and operators . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.4.1 Logical operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.4.2 Relational operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.4.3 Logical indexing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.5 Cell arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

iv
 3.6 Structure arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.6.1 Dot-notation and nested structures . . . . . . . . . . . . . . . . . . . . 16
3.6.2 Indexing structure arrays . . . . . . . . . . . . . . . . . . . . . . . . . 16

4 Structured programming 17
4.1 Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.1.1 Writing your own functions . . . . . . . . . . . . . . . . . . . . . . . . 17
4.1.2 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2 Conditional statements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2.1 if, elseif, else . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2.2 switch, case, otherwise . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2.3 try, catch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.3 Loops and iterations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.3.1 while . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.3.2 for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.4 Control statements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5 Basic graphics functions 21

5.1 The plot command . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2 Customising the Axes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2.1 Axes scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2.2 Closing Figure windows . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2.3 Adding lines to existing Axes . . . . . . . . . . . . . . . . . . . . . . . 22
5.2.4 Axes labels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.2.5 Adding a Legend . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
5.3 Exporting to file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.3.1 Raster vs. vector graphics . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.4 Creating additonal Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5.5 The subplot command . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5.6 Other graphics commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

6 Advanced graphics programming 26

6.1 Graphics objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.1.1 Top-level objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
6.1.2 Child objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
6.2 Object handles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
6.2.1 Get current . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
6.3 Object properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
6.3.1 Setters and getters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
6.3.2 Dot notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
6.3.3 Defining properties at declaration . . . . . . . . . . . . . . . . . . . . 29
6.4 Creating a publication-quality Figure . . . . . . . . . . . . . . . . . . . . . . . 29

7 The Control System Toolbox 32

7.1 Dynamic system models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
7.1.1 Transfer function models . . . . . . . . . . . . . . . . . . . . . . . . . 33

v
 7.1.2 Zero-pole-gain models . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
7.1.3 State-space models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
7.1.4 PID controller models . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
7.2 Interconnecting dynamic systems . . . . . . . . . . . . . . . . . . . . . . . . . 36
7.2.1 Parallel connection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
7.2.2 Series connection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
7.2.3 Feedback connection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
7.3 Linear analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
7.3.1 Time-domain analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
7.3.2 Frequency-domain analysis . . . . . . . . . . . . . . . . . . . . . . . . 39
7.3.3 Pole-zero analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

8 Dynamic systems – Advanced topics 43

8.1 Model objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
8.2 Multi-variable systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
8.2.1 State-space models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
8.2.2 Transfer matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
8.3 Discrete-time systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
8.3.1 The discrete-time domain . . . . . . . . . . . . . . . . . . . . . . . . . 46
8.3.2 Declaring discrete system models . . . . . . . . . . . . . . . . . . . . 47
8.3.3 Converting discrete system models . . . . . . . . . . . . . . . . . . . . 47
Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
8.4 Customisable plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

9 Simulink 50
9.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
9.2 Fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
9.2.1 Getting started . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
9.2.2 Library Browser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
9.2.3 Running a simple example . . . . . . . . . . . . . . . . . . . . . . . . 53
9.3 Solver configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
9.3.1 Fixed-step methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
9.3.2 Variable-step methods . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
9.3.3 Stiff systems of equations . . . . . . . . . . . . . . . . . . . . . . . . . 55
9.3.4 Discrete solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
9.3.5 Automatic selection by Simulink . . . . . . . . . . . . . . . . . . . . . 56
9.4 Studying a demo model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
9.4.1 Coordinate transforms . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
9.4.2 Subsystems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
9.4.3 Simulation resuls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
9.5 Creating your first dynamic model . . . . . . . . . . . . . . . . . . . . . . . . 61
9.5.1 Walkthrough . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
9.5.2 Partitioning the Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

vi
 9.6 Matlab integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
9.6.1 Exporting outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
9.6.2 To Workspace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
9.6.3 Running simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
9.6.4 Single simulation output . . . . . . . . . . . . . . . . . . . . . . . . . . 69
Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
9.7 Dynamic system models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
9.7.1 The Transfer Fcn block . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
9.8 Closing words . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Solutions to exercises 72
Elementary mathematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
Vectors and matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
The Control System Toolbox . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Discrete-time systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
Simulink-Matlab integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

Appendix 80
Creating a publication-quality Figure . . . . . . . . . . . . . . . . . . . . . . . . . . 80
Custom linear analysis plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

vii
1 Introduction

Matlab is one of the most widely used languages for scientific programming and numerical
computation, with over 3 million users worldwide. Typical applications include, but are not
limited to data analytics, deep learning, signal processing, robotics, and control systems.
It is an interactive language, meaning that it will execute commands directly, without
first being compiled into machine code. The fundamental datatypes are double-precision
floating point matrices (Matlab stands for MATrix LABoratory). This is one of the reasons
that it performs so well in the domain of numerical computing.
“Matlab” also refers to the application itself, which is required to run code written in the
language. It is an Integrated Development Environment (IDE) that can handle all necessary
aspects of developing sophisticated applications in the Matlab language. In addition to
the interactive command line, it provides tools for graphical drawing and writing your own
scripts and functions, to name a few features. This tutorial will focus on both the language,
regarding syntax and semantics, and the IDE itself.
Simulink is a graphical user interface that sits on top of the Matlab engine and simplifies the
modeling and simulation of complex dynamical systems. The necessary matrix calculations
are performed by Matlab in the background, automatically. Both products are proprietary
(non-free) and are developed and maintained by the US-based company MathWorks, Inc.
Please note: The tutorial is written with MATLAB R2021a in mind. In the case that you
are working with another version, some functions may behave differently, and some menus
and dialogues might have a different appearance.
In addition to the core packages Matlab and Simulink, you will need to have the following
packages installed for this tutorial
• Control System Toolbox (Chapters 7+)
• System Identification Toolbox (Chapter 8)
When installing Matlab, it is suggested that you include all packages, if you have the
capacity.

1
2 Fundamentals

2.1 Introducing the desktop

After starting the Matlab application you will be greeted with the desktop user interface,
see Figure 2.1. It consists of multiple panels, where the Command Window with the
command line “>>” dominates the screen. A brief description of the panels is given below.

Figure 2.1: Default appearance of the Matlab user interface

• Current Folder – Access the files and folders in the current working directory
• Command Window – Enter commands at the command line, more on this later
• Workspace – Explore variables and data that you create or import

Next, consider the elements located above the panels.

• Toolstrip – Contains quick access to creating or opening scripts, changing settings, etc.
• Navigation Toolbar – Allows changing the current working directory

2
 2 Fundamentals

2.2 The command line

Now to get started. The Command Window serves as the standard input-output between
the user and the Matlab engine. Every command provided by the user is completed with
Enter , upon which any output will be returned below it.
Matlab is an inferred-type language, meaning variables can be created without explicitly
declaring their type (unlike C or Java for example). Try entering the following commands
into the Command Window and see if you understand what is happening
>> 6*7
>> a = 3* pi
>> b = a/5
>> c = a*b
>> b = 1:3:20
>> plot(b)
>> help cos
>> b = 2cos(b)

The last line will trigger one of the most common error messages:
Error: Invalid expression ...

Matlab will sometimes try to help you find the cause of your problems. Try correcting the
command according to the advice you get, and see if you eliminate the error.

2.2.1 Suppressing output

The output in the command window can be suppressed with a semi-colon at the end of the
line. Try this out by repeating the above commands with an added semi-colon. But first,
read on!

2.2.2 Input history

Instead of typing the commands anew, you can use the Up Arrow to navigate through
the most recently used commands. Any previously entered line can be changed before
being entered again with . You can find the previous command even faster if you type
the first few letters of the line before pressing – only matching lines will be shown to
you then.

2.2.3 Auto-completion
Finally, to speed things up in the future, it is convenient to use Tab to auto-complete
your typed commands. Test this by creating a variable that shares your name, and is equal
to your age (in years). Now start typing the first 2-3 letters of your name in the command
line, then press to complete the rest. Depending on your name, it may prove more
or less useful, but for structured variable names such as motorModelLinearSimplified it
may help. Note that it will look for completions among the default functions and variable
names, including other toolboxes that are installed on your system.

3
 2 Fundamentals

2.2.4 Documentation
Oftentimes, you will need to look up a function to see exactly what syntax it uses, which
input arguments are needed, or what to expect it to return. To facilitate this, the documen-
tation for functions can be accessed through the command line directly. Try the following
commands
>> help cos
>> help linspace
>> doc linspace
>> lookfor cosine

The command lookfor searches through the first line of all documentation contained in
the search path, to try to find any and all matches for the query. Searching through the
documentation may sometimes take longer than you wish for, and the Matlab window
will be unresponsive while this happens.

2.2.5 Aborting commands

To abort any running command, you can use Ctrl + C .

2.2.6 Unquoted strings

As you will have noticed, the second word that you type is highlighted in purple. This is
not just an aesthetic detail, but rather, is indicative of what is happening. Try running the
following commands
>> disp Test
>> disp('Test ')

Then try to run the following

>> disp Hello World
>> disp 'Hello World '

Based on the error message you receive, you might be able to deduce what is happening.
If you write unquoted strings in Matlab, it will parse the first word as the function, and
any additional words separated by spaces as separate input arguments. By using single
apostrophes, you can forcefully tie together strings of letters.
It should be said that when writing scripts, it is considered good practice to explicitly use
brackets when calling functions to improve readability.

2.2.7 Housekeeping
You can tidy up in Matlab using the following commands
>> clear all % Clear all variables from the Workspace
>> clc % Clear the Command Window

4
 2 Fundamentals

2.3 Scripts
When you wish for several commands to be performed successively, you can write your
Script, also sometimes called M-Files because of their “.m” ending.
When running a script, all variables defined in the (base) Workspace will be available to
it, and conversely, all new data created by your script will be stored there as well. These
variables are said to be in the global scope. This is a significant difference when compared
to Matlab functions, which are addressed in Section 4.1. Functions, unlike scripts, each
have their own scope and are not able to access the Workspace directly.

2.3.1 Creating new scripts

To create a new script, click the “New Script” button located in the top left of the Toolstrip,
or type
>> edit

This will open the Editor window, which will appear above the Command Window. The
Editor has syntax highlighting, supports debugging, and enables formatting code using
the right-click context menu.
Observe that the Toolstrip now has a new tab, labeled “EDITOR”, providing additional
tools for code formatting. Once you have populated your script with the command you
want to run, you can click “Run” in the Toolstrip, or F5 . This will prompt you to save your
script somewhere, and the file is saved every time you run it this way.
Create a new script, and write the following in the Editor
% My first script
cd
disp('This is my first script ')
a = 6*7
b = 6*7;

Save it as “myScript1.m”.
Then try the following commands in the command line
>> myScript1
>> help myScript1

2.3.2 Sectioning
For long scripts, it can sometimes be efficient to divide the code into Sections, which helps
break down the code both visually but also enables running parts of it at a time.
This is done by inserting Cell Dividers (%%) above the parts you want to section. The section
where the line marker is currently in will also be highlighted in yellow. Sectioning allows
for quick traversal through your code by using the keyboard combinations Ctrl + and
Ctrl + .

5
 2 Fundamentals

You can run the current section (without saving) with Ctrl + . This will also ignore any
Breakpoints in the code, which will be explained in the next section.

2.3.3 Debugging
When you want to assure yourself of what is happening in the code, you can place Break-
points in the code, by clicking the line number in the Editor. A red circle will show itself,
which represents a breakpoint in the code. Every time the execution encounters this line,
it will pause there, and the Workspace will reflect the state it was in at that time.
When the script has stopped for a breakpoint, you will see the message “Paused in debug-
ger” in the bottom left, and the command line will look like
K>>

You can use the command line during this state to manipulate the Workspace, to test the
behavior in the code. To stop debugging, click “Quit Debugging” in the Toolstrip, or use
+ F5 .

6
3 Basic functions and datatypes

This chapter aims to introduce some of the building blocks of writing code in Matlab, and
we will begin by introducing some of the built-in functions for mathematical operations.
After this, you will be introduced to the two most important datatypes in Matlab – the
numeric and the logical, and the notion that everything is a matrix. Finally, we will look at
the structure array, which allows you to create your own data structures.

3.1 System variables

As shown in the previous chapter, you are free to declare as many variables as you like, but
it should be said that Matlab contains some pre-defined, system variables (see Table 3.1)
whose names you should refrain from using.

Table 3.1: Some of the reserved system variables in Matlab

ans Contains the result returned from the most recent command
pi An approximation of the number 𝜋
eps Machine epsilon, floating-point relative accuracy, 2−52
Inf Infinity, represents 1/(+0)
NaN Not-a-Number, represents 0/0
√
i, j The imaginary unit, −1

3.2 Elementary mathematics

Matlab is developed especially for mathematical and scientific programming, and as such,
contains a rich catalog of mathematical functions and operations that are pre-defined in
every installation.
Unlike other programming languages, no library needs to be “included”, or “imported”,
all built-in functions are available from the start. To see a list of some of the elementary
mathematics functions, you can type
>> help elfun

An excerpt of this list is shown in Table 3.2.

7
 3 Basic functions and datatypes

Table 3.2: Some of the elementary mathematics functions in Matlab

sin Sine sqrt Square root

cos Cosine abs Absolute value
exp Exponential mod Modulus
log Natural logarithm abs Signum

Exercise
Create a new script, and save it as EulersFormula.m and write the following code

1. Declare the variable 𝑥 = 𝜋

2. Write the following two expressions in Matlab, using elementary functions and
variables

𝑓1 = cos 𝑥 + 𝑖 sin 𝑥 𝑓2 = 𝑒 𝑖 𝑥

3. Show that the answers of the two expressions are equal for 𝑥 = 𝜋, by showing that
the difference is zero
4. Try some other values of 𝑥, run the script, and see what happens

Using the same script, investigate the following

1. Let 𝑥 = 1
2. Write the following two expressions
 
1 1 1 1
𝑔1 = 𝑔2 = √ + √
1+𝑥 2 1+𝑖 𝑥 1−𝑖 𝑥

3. Show that the answers of the two expressions are equal for 𝑥 = 1
4. Try the following values: 0, −1, 2, 𝑖 + 2 and try to explain what you are seeing

8
 3 Basic functions and datatypes

3.3 Vectors and matrices

Matlab is developed for use in matrix calculations. A matrix is considered by Matlab
to be a field of numerical entries, represented by double-precision floating-point decimals.
Strings, such as 'this' for example, is just a vector containing characters1 . Scalars are
represented internally as 1 × 1 matrices, vectors as 1 × 𝑛 matrices (column vectors), or 𝑛 × 1
matrices (row vectors).
When using matrix operations (like multiplication), you have to take into consideration
what you want to calculate, and what kind of matrices you have to use, otherwise, you
will likely get error messages, or worse, end up calculating something else than what you
wanted.

3.3.1 Declaring matrices

Matrices are declared using square brackets [ ], where each entry is separated by a comma
and each new row by a semicolon, like this
>> u = [1 ,2 ,3] % Row vector
>> v = [4;5;6] % Column vector
>> A = [1 ,2 ,3;4 ,5 ,6;7 ,8 ,0] % Matrix
>> c = 'test ' % Character vector

Row vectors with fixed intervals can be created using : (colon)

>> k = 1:1:10 % Initial value : Step -size : Final value

Alternatively, you can create a vector with a linear distribution using

>> x0 = 0; % Initial value
>> x1 = 2* pi; % Final value
>> n = 100; % Number of elements
>> x = linspace (x0 ,x1 ,n)

Other special matrices can be created using the following commands

>> u = zeros (1 ,5) % Matrix populated with only zeros
>> v = ones (2 ,4) % Matrix populated with only ones

To see a more comprehensive list of elementary functions for declaring and manipulating
matrices, type
help elmat

1 Starting in R2017a, you can create strings using double quotes, which provide a set of functions for working
with text as data. Read more on the Mathworks homepage

9
 3 Basic functions and datatypes

3.3.2 Diagonal matrices

A diagonal matrix has non-zero entries only along its diagonals. The notion is extended to
include upper and lower diagonals, such as the bi- or tridiagonal matrix, and so on.
The most common diagonal matrix is the identity matrix, which is a square matrix with
only ones along the main diagonal.
>> I = eye (3) % The identity matrix

You can create square diagonal matrices using diag.

>> d = 1:1:3; % Declare 1x3 vector
>> D = diag(d) % Create 3x3 diagonal matrix

You can specify which diagonal using an extra integer argument. Giving no argument is
implicitly zero, and refers to the main diagonal.
>> g0 = 2* ones (1 ,4); % The main diagonal
>> g1 = 3* ones (1 ,3); % The first upper diagonal
>> G = diag(g0 ,0) + diag(g1 ,1);

3.3.3 Sparse matrices

A more versatile way of declaring diagonal matrices with several diagonals is to use the
spdiags. The “sp” is for sparse, which is a lightweight way of storing large matrices with
few non-zero elements in memory. The interesting aspect here is not only the sparse nature
of the resulting matrix but the flexibility of the command syntax.
>> n = 5; % Number of rows and columns
>> h = ones(n ,2); % Diagonals grouped in columns
>> k = [0 ,1]; % Main and first upper diagonal
>> H = spdiags (h,k,n,n); % Create sparse matrix
>> full(H) % Display full storage form

It should be mentioned that you can freely convert between sparse and full representa-
tions; numerically there is no difference between the two.

3.3.4 Accessing elements

If you want to access individual elements in matrices you can use the desired row and
column index, inside regular brackets ( ). The indices themselves can also be vectors, to
access several rows or columns at once.
>> A = [1 ,2 ,3; 4 ,5 ,6; 7 ,8 ,9]
>> A(1 ,3) % Element in row 1, column 3
>> A(1 ,:) % ':' here represents the entire column
>> A(: ,1) % ':' here represents the entire row
>> A([1 ,3] ,1) % Elements 1 and 3 from column 1

10
 3 Basic functions and datatypes

Note that Matlab uses one-based indexing, unlike some languages that start counting at
zero. For vectors, you can drop the second index, which is then implicitly 1. You can also
use the reserved keyword end to access the last entry in a matrix.
>> b = 1:1:6
>> b(end -1) % Second -to -last element
>> b(4: end) % Last three elements
>> b(end +1) = 7 % You can also add elements like this
>> b(3) = [] % And remove elements with an empty matrix

You can also return the dimensions of declared matrices using the functions size and
length, try the following commands
>> size(A)
>> length (A)
>> size(b)
>> length (b)

3.3.5 Matrix operations

Most of the time, operators defined in Matlab are defined for matrices and will try to inter-
pret the result as a matrix operation. This includes things such as division, multiplication,
power, etc. All matrix operations must obey the mathematical rules governing matrices;
for instance, matrix multiplication of two matrices, 𝐴 and 𝐵, is only defined if the number
of columns of 𝐴 is equal to the number of rows of 𝐵. Similarly, you cannot add together
two matrices of different dimensions.
Lastly, matrix multiplication and “division” is generally not commutative, as is the case
with scalars. Matrix division in Matlab is performed by solving the linear system of
equations, and the equivalent mathematical expressions for invertible matrices are shown
below

𝐵/𝐴 = 𝐵 𝐴−1 𝐴\𝐵 = 𝐴−1 𝐵

Table 3.3 shows a list of some of the commonly used arithmetic operations for matrices and
their corresponding operator.

Table 3.3: Some of the matrix operators in Matlab

mtimes * Multiplication
mpower ^ Power
mldivide \ Left division, backslash
mrdivide / Right division, slash

11
 3 Basic functions and datatypes

3.3.6 Element-wise operations

Sometimes you find yourself not wanting to perform matrix operations at all, but you want
to apply an operation element-wise to each element in a matrix. Table 3.4 shows a list of some
of the commonly used element-wise operations. Note that some of them are variations of
the matrix operations mentioned earlier.

Table 3.4: Some of the element-wise operators in Matlab

plus + Addition
minus - Subtraction
times .* Element-wise multiplication
power .^ Element-wise power
ldivide .\ Element-wise left division
rdivide ./ Element-wise right division
transpose .’ Transpose

To get a feeling for the difference between the two kinds of operators, try the following
>> A = [1 ,2 ,3; 4 ,5 ,6; 7 ,8 ,9] % Declare a 3x3 matrix
>> b = [5 ,6 ,7] % Declare a 1x3 vector
>> b = b.' % Transpose to a 3x1 vector
>> b + b % Vector addition
>> b.' * b % Scalar (dot) product
>> b * b % Error : Incorrect dimensions !
>> b .* b % Element -wise multiplication
>> A * A % Matrix multiplication
>> A ^ 2 % This gives the same result
>> A .^ 2 % But not this!
>> A \ b % Backslash division

You can see a comprehensive list of all the arithmetic operations, and more, if you type
>> help ops

3.3.7 Matrix functions

Finally, there is an extensive library of built-in functions for numerical linear algebra and
matrix analysis, which can be seen with the command
>> help matfun

Some examples are shown in Table 3.5

12
 3 Basic functions and datatypes

Table 3.5: Some of the matrix analysis functions in Matlab

rank Matrix rank eig Eigenvalues and eigenvectors

det Determinant svd Singular value decomposition
inv Matrix inverse expm Matrix (not element-wise) exponential

Exercise
Very rarely will you declare matrices by hand; oftentimes, you will want to build your
matrices dynamically, varying only a few parameters to shape the outcome. To wrap up
this section, you will try your hand at creating an automatic build script for a generic matrix.
For this exercise, you will create the following sparse, diagonal, 𝑚 × 𝑛 rectangular matrix.

 𝑎0 𝑎1 · · · 𝑎𝑘 
𝑎0 𝑎1 · · · 𝑎𝑘
 
𝐴 = 
 
.. .. .. .. 
 . . . . 


 𝑎0 𝑎1 ··· 𝑎 𝑘 

The goal is to automate all the work, without using any “hard-coded” constants. The user
need only provide the vector to be repeated along the diagonal, and the number of rows.
Do the following in a new Matlab script.
1. Declare the vector 𝑎 = [1, 2, 3] and choose 𝑚 = 3.
2. Find a way to automatically determine 𝑛, the number of columns of the matrix 𝐴.
3. Using 𝑎 and 𝑚, create a matrix whose columns represent each non-zero diagonal.
Hint: To repeat a matrix several times and concatenate the result, check out
>> help repmat

4. Next, declare the appropriately sized index vector for the diagonals.
5. Finally, create the sparse diagonal matrix 𝐴 using spdiags. You should convert it to
its full form to more easily see if it looks correct.
6. Now try to change 𝑚 to some other positive integer and run the script again. If you
have done it correctly, the script should produce a new, appropriately sized matrix –
without any complaints.
7. Then, change the vector 𝑎 to some other, non-empty row vector of any size. Again, if
the process is automated properly, there should be no error.
8. The “trivial” case where 𝑎 = 𝑚 = 1 should also work.

13
 3 Basic functions and datatypes

3.4 Logical values and operators

It is sometimes interesting to compute not the numerical, but the logical value of an expres-
sion. In other words, determining whether an expression is true or false.
In Matlab, the datatype that stores logical values is simply referred to as a “logical”, and
like numerics, can be a scalar, vector, or matrix. The logical can only assume one of two
values; logical 1 (true), or logical 0 (false). Logical matrices can be directly declared using
the functions true and false, but you can also cast numbers into logicals, and it is defined
such that any non-zero number is a logical 1, and only zero is a logical 0. Try the following
>> x = -1:1:2
>> logical (x)
>> true
>> false (3)
>> true (1 ,5)

3.4.1 Logical operators

It is possible to test the logical value of matrices using logical operators, such as AND, OR,
and NOT, shown in Table 3.6.

Table 3.6: Some of the logical operators in Matlab

and & Element-wise AND

or | Element-wise OR
not ~ Element-wise NOT

For example
>> A = [5 ,7 ,0; 0 ,2 ,9; 5 ,0 ,0]
>> B = logical ([1 ,1 ,0; 1 ,0 ,1; 1 ,0 ,0])
>> A & B % All true or non -zero elements in A AND B
>> A | B % All true or non -zero elements in A OR B
>> ~B % All false elements in B
>> ~(A | ~B) % you can also write compound expressions

3.4.2 Relational operators

To compare the numerical values of two matrices, you can use the relational operators, a
few are shown in Table 3.7.
For example
>> B = [6 ,6 ,0; 1 ,3 ,5; -1,0,0]
>> A == B % All elements equal in A and B
>> A > B % All elements in A greater than B
>> B <= 3 % All in B less than or equal to 3

14
 3 Basic functions and datatypes

Table 3.7: Some of the relational operators in Matlab

eq == Determine equality
ge >= Determine greater than or equal to
gt > Determine greater than
le <= Determine less than or equal to
lt < Determine less than
ne ~= Determine inequality

3.4.3 Logical indexing

It is also possible to use local matrices when indexing matrices. The result is a vector
containing all elements whose index corresponds to a logical 1, removing all other values.
Similarly, you can multiply any two valid numeric and logical matrices, in which case the
logical is cast into a numeric one or zero.
>> X = (A >= 2) & (A <= 5) % All elements a, s.t. 2 <= a <= 5
>> A(X) % Extract the elements in a vector
>> A .* X % Set all false values to zero

3.5 Cell arrays

In Matlab a Cell array is a generalized class of indexed data containers, called cells. Each
cell can contain any type of data, such as text, matrices of different sizes, or even other cell
arrays.
Cell arrays are declared using curly braces { }, as shown below.
>> C = {1, 2, 3; 'text ', [3 ,4 ,5] , {6 ,7 ,8} }

You can access subsets of cells using regular brackets ( ).

>> C(2 ,1) % Cell located in row 2, column 1
>> C(1 ,1:3) % Subset of cells located in row 1

You can access the contents of cells by indexing using braces.

>> C{2 ,1} % Element contained in row 2, column 1
>> C{1 ,1:3} % All elements contained in row 1

Note that when accessing multiple elements contained inside cells, the output is returned
individually, not as a matrix. Cells are often used to store text strings, since storing them in
the same matrix would require vertical concatenation (stacking them on top of each other),
and so they would need to have the same length.

15
 3 Basic functions and datatypes

3.6 Structure arrays

To create more complex datatypes, you can use structure arrays, or structs for short. A
struct is a generic container that can be filled with an arbitrary amount of fields, which can
contain data of different types. To give an example, consider that we are simulating some
point 𝑝 1 in a 3D space, which is defined by its coordinates [𝑥, 𝑦, 𝑧].
Let’s represent the coordinates as regular 1×3 vectors, and assume that the initial conditions
are zero.
We can then create the appropriate struct using
>> p = struct ('pos ', zeros (1 ,3))

Now assume we also want to keep track of the point’s velocity, [𝑥,
¤ 𝑦,
¤ 𝑧¤ ]. Let’s give it an
impulse along the 𝑧-axis. One way to add this new information is to use the setfield
command, as follows
>> p = setfield (p, 'vel ', [0 ,0 ,1])

Now if you double-click the variable p in the Workspace, it will bring up the Variables
window, and show the fields contained inside it.

3.6.1 Dot-notation and nested structures

Next, let’s give the point a header, which will help us keep track of its name and anything
else which could be important. We can store this inside a new struct, contained in the
original struct.
Let’s also introduce a more efficient way of accessing fields, using . (dot) notation.
>> p. header = struct ('name ', 'Sputnik ', 'id ', 2)

As can be seen, the dot-notation can be used even if the field does not yet exist, upon which
it will be created. It can also be used multiple times to go deeper into the data structure;
for instance, you can access the name of the point using
>> p. header .name

Usually, a header comes first, so let’s sort the order of the fields using
>> p = orderfields (p) % No extra arguments -> sort alphabetically

3.6.2 Indexing structure arrays

Now say that we want to add one more point. We could do the same thing again, and
create a new struct containing the same fields. Or, we can use the fact that structs can be
indexed, and simply add another entry in the structure array 𝑝:
>> p(2) = struct ('vel ', [1,-1,0], 'pos ', [1 ,1 ,0] , 'header ',
struct ('name ', 'Cassini ', 'id ', 25008) )

Note that it does not matter in which order fields are specified.

16
4 Structured programming

This chapter covers some advanced techniques and constructs that, while they are useful
when developing code, lie outside the scope of modeling mechatronic systems. The reason
this chapter is included is more for the sake of completeness.

4.1 Functions
One way to help structure your program is to separate the project code into functions, which
also allows for reusing functionality instead of re-writing it.
Another benefit is that the code becomes easier to document for the developer and easier
to read for an external user. Note that this could be you in a few years, re-visiting an old
project.

4.1.1 Writing your own functions

A function works much like a script, and you create it in the same way that you would a
normal M-file, but it has some key differences. Consider the following example definition
of the function myfunc.
function <output> = myfunc ( <input> )
% MYFUNC A brief summary of the function .
% MYFUNC (X) An example use -case , what to expect from X.
%
% Additional information .

...

return
end

First, note the signature in the first line of code, which lets Matlab know to treat the file
as a function definition, and how to handle the ingoing and outgoing values. Next, the
function definition makes use of the reserved keyword return, which stops the execution
of the function and returns control to the place where it was invoked. Finally, the block of
code is enclosed by the keyword end, which marks the end of the function definition.
Create a new (empty) M-file, and add the following lines of code.

17
 4 Structured programming

function result = add(a, b)

%ADD Addition of two numbers .
% ADD(X,Y) returns the sum of the numbers X and Y.
%
% This is my first function .
result = a + b;
return
end

Save it as add.m. Observe that it will appear in the Current Window with a slightly different
icon (“fx”) to let you know it’s a function definition.
Now try the following commands
>> help add
>> add (2 ,2)

4.1.2 Scope
Note that in the previous example, no variables called a or b ever appeared in the base
workspace. An important aspect of functions is that they have their own private workspace,
which is destroyed at the end of the function call. This also means that a function cannot
directly access the base workspace, unless specified explicitly using the global keyword,
both inside the function definition, and in the scope of the invoking function.

4.2 Conditional statements

Like in all other programming languages, Matlab enables the use of logical conditions
to selectively alter the control flow of the program. This section presents examples of the
Matlab syntax.

4.2.1 if, elseif, else

The if-statement is the most fundamental of all the conditional statements. It makes the
program execute a subset of the defined commands, based on the logical value of some
condition, which is placed inside the corresponding set of brackets.
It is possible to have several conditions available; if one of the conditions is satisfied, the
program will execute the specified subset of commands, then go immediately to the line
after the end keyword, ignoring the other options.
if <condition>
% Go here if the first condition is satisfied
elseif <condition>
% Go here if only the second condition is satisfied
else
% Go here if none of the above conditions are satisfied
end

18
 4 Structured programming

4.2.2 switch, case, otherwise

The switch-statement is used when you want to control the program based on the value
of a single variable or expression. This is just a special case of the above statements, where
each condition tests different values of the same expression.
switch <expression>
case <value>
% Go here if the expression equals the first value
case <value>
% Go here if the expression equals the second value
otherwise
% Go here if the expression equals something else
end

4.2.3 try, catch

Sometimes, exceptions can arise during the execution of your code, which can cause un-
expected behavior. Instead of crashing the execution and producing an error, it is often
preferred to provide customized feedback to the user, close any open data streams, and
terminate the program gracefully. This can be achieved with the try-catch-statements.
Commonly, you would put commands that are known to be sensitive inside such a state-
ment; this is often used where a mistake in user input may cause undefined behavior, or
where the program depends on some resource that might not always be available. The
exception can be stored in the workspace by specifying a variable after the catch-keyword
try
% Try something that normally works
catch <exception>
% Do something special in the case of an error
end

It is sometimes possible to completely recover from exceptions and continue execution, by

redirecting the control flow (see Section 4.4).

4.3 Loops and iterations

Similar to conditional statements, a loop is a control structure that can redirect the flow
of the execution, based on certain conditions during runtime. Just like the conditional
statements, the block of code consists of everything between the loop statement and the
keyword end.
Any statements defined inside the block can be executed several times, once or never, based
on the type of function used.

19
 4 Structured programming

4.3.1 while
In a while-loop, the specified statements will be continuously executed until the condition
resolves to a logical zero. In the case where the condition resolves to false the first time, the
loop will not be entered, and the statements will never be executed.
while <condition>
% Do something repeatedly , while the condition holds
end

If the logical value of the condition is not altered by any commands inside the loop, the
program can get stuck in an infinite loop. Such a loop can be created intentionally with the
condition while true, which will never be exited unless a control statement is used, see
Section 4.4

4.3.2 for
The for-statement is suited for cases where the number of iterations is defined beforehand.
The index variable will sequentially assume each defined value, and the code will be
executed once per iteration. When all values have been used, the program will break out
of the loop, and keep executing.
for <index> = <values>
% Do something iteratively , for each specified value
end

4.4 Control statements

It is possible to manipulate the program execution further by using the special, reserved
keywords listed below.

• break can be used to immediately break the program out of the current loop and
execute commands defined after the end of the loop.
• continue is used to immediately go to the next iteration of the loop.
• return can be used to exit the current function, and return control to the invoking
function. If this is used outside of a function it will return control to the command
prompt, causing the program to terminate.

20
5 Basic graphics functions

One of the things that makes Matlab so useful is its rich library for graphical representation
of data, which helps to visualize data and communicate results. The built-in graphics
functions include plotting functions in 2D and 3D, animating objects, applying lighting or
transparency effects, and adding user interface elements and callbacks.
This chapter will show how to create graphical representations of data interactively, using
higher-abstraction functions that simplify some of the features behind the scenes. It will
cover some 2D plotting functions, modifying and working with multiple figures, and
exporting to vector graphics formats. If a larger amount of control is desired, it is possible
to use the lower-level functions to build up graphical objects from the start, which will be
elaborated on in Chapter 6.

5.1 The plot command

As an introduction, let’s look at the most commonly used graphics command; namely plot.
Try the following commands
>> x = linspace (0 ,2*pi ,30);
>> y1 = sin(x);
>> y2 = cos(x);
>> plot(y1) % Compare this ...
>> plot(x, y1) % with this
>> plot(x, y1 , y2) % This won 't work ...
>> plot(x, y1 , x, y2) % but this will!

As you can see, the default behavior of the plot command is to open a new figure, populate
it with axes, scale them accordingly, and draw lines onto it, making the first line blue, then
red, etc. Sometimes this is sufficient, but it is also possible to control the appearance, as
will be shown in the coming sections.

5.2 Customising the Axes

The word “axis” refers to one of the dimensions of a plot; the 𝑥-axis for instance, and the
plural form is “axes”. The word “Axes” (note the capitalization) in this text refers to the
Axes object, which is the canvas upon which data is visualized. This object contains the
(𝑥, 𝑦)-axes, including the tick marks, text labels, grid lines, and the background color. You

21
 5 Basic graphics functions

can set some basic Axes properties using the keywords 'on' or 'off', try the following
commands
>> grid on % Show grid lines
>> box on % Enclose the Axes in ruler lines

5.2.1 Axes scaling

You can configure the Axes size and scaling using the axis command. Look at the docu-
mentation (help axis) and try some of the possible commands, such as
>> axis % Return the current scaling as a vector
>> axis ([3 , 5, -1, 0]) % Set custom scaling
>> axis tight % Set the axis limits to the range of the data
>> axis equal % Set the aspect ratio for tick marks to 1:1
>> axis off % Hide the Axes

5.2.2 Closing Figure windows

You can close open Figures with the following command
>> close % Close the current Figure
>> close all % Close all open Figures

5.2.3 Adding lines to existing Axes

A call to plot will remove existing lines from the current Axes. To disable this you can set
the hold property of the Axes.
>> plot(x, y1) % Create a new plot
>> hold on % Hold on to the current lines
>> plot(x, y2) % Add an additional line
>> hold off % Disable hold

5.2.4 Axes labels

To add text labels to your Axes, use the following commands
>> xlabel ('Angle (rad)') % x-axis label
>> ylabel ('Amplitude ') % y-axis label
>> title ('Trigonometric functions ') % Axes title

5.2.5 Adding a Legend

In the case where you have multiple lines in the same plot, it must be clear what you are
visualizing. To provide sufficient information to interpret the data, it is often necessary to
add a Legend to your Figure.
>> legend ('Sine ','Cosine ')

The strings should be provided in the order that the lines were added.

22
 5 Basic graphics functions

5.3 Exporting to file

The most user-friendly way to export your figure to a document is to use the File Save As
dialog in the Figure window. You can then choose to save the Figure to various formats,
the most common are shown in Table 5.1.

Table 5.1: Some of the supported export formats supported in Matlab

fig Matlab Figure Internal format

bmp Microsoft Bitmap Raster graphics
jpg JPEG Image Raster graphics
png Portable Network Graphics Raster graphics
tif Adobe Tagged Image File Format Raster graphics
eps Encapsulated PostScript Vector graphics
pdf Portable Document Format (full page) Vector graphics
svg Scalable Vector Graphics File Vector graphics

The .fig format is used for saving the Figure as-is, so that it can be loaded at some later
stage and further manipulated inside Matlab. The other formats involve an element of
rendering, in which an image is produced, which is meant to capture what is displayed on
the screen. The image formats can be divided into two categories; namely raster (also called
bitmap), or vector graphics.

5.3.1 Raster vs. vector graphics

The fundamental difference between the two formats is that raster images are composed
of a grid of pixels, and so are resolution dependent. This means that if you zoom into
a raster image beyond a certain point, you will start seeing the pixels that compose the
image. Some raster images also include destructive compression algorithms, which can
irreversibly degrade the image quality.
An image in a vector graphics format, on the other hand, is rendered by the document
viewer for the current screen resolution. When you change the zoom level, the view of the
document is rendered again to produce a new view of the document. Regardless of how
far you zoom in, the graphics or text is not going to become pixelated.
Compare the two images in Figure 5.1. They were both exported through the “Save As”
dialogue in Matlab, using two different image formats. Try different zoom levels when
comparing the images, and try to highlight the text. Which one looks the best?
The vector graphics file in the example is exported to the Encapsulated PostScript (EPS)
format, which is a widely used vector graphics format, supported by the graphicx package
in most modern LATEX distributions, like MiKTeX, TeX Live, and Overleaf.
The raster image is rendered using the JPEG format, which implements a lossy compression
algorithm, which adds to the reduction in image quality. While other raster formats may

23
 5 Basic graphics functions

Trigonometric functions
1

0.8

0.6

0.4

0.2

Amplitude
0

-0.2

-0.4

-0.6

-0.8

0 1 2 3 4 5 6
Angle (rad)

(a) JPEG: Compressed bitmap (20 KB) (b) EPS: Vector graphics (18 KB)

Figure 5.1: Comparison between a compressed bitmap and a vector graphics print

provide less destructive compression, it helps to illustrate one of the potential issues. It
should be said that there are a few use cases where raster graphics is preferred, most
commonly when plotting 3D surfaces, or colors with transparency.
If your vector graphics file still appears pixellated, make sure you are using the painters
renderer. You can configure this under “Export Setup...” in the Figure window. Go to the
“Renderer” tab and change the “Custom renderer” option.

5.4 Creating additonal Figures

Up until now, any plotting commands have been placed in the already open Figure, with
one being created only if there is not one open. You can create additional Figures with the
figure command. Any subsequent commands will target the most recent window
>> figure % Open a new Figure window
>> plot(x, x.^2) % Populate it with a new plot

5.5 The subplot command

Sometimes it is desirable to display multiple Axes inside the same Figure, when the plotted
quantities have different units or are orders of magnitude apart, for example.
You can divide the current figure into a 𝑛 × 𝑚-grid by using subplot(n,m,p), where any
subsequent plotting commands are applied to the Axes in the position specified by 𝑝.
>> figure % Open a new Figure window
>> subplot (2 ,2 ,1); % Create a 2x2 Axes grid
>> plot(x, y1); % Plot in the first grid element
>> title ('Subplot 1: Sine ')
>> subplot (2 ,2 ,2); % Select the second grid element
>> plot(x, y2);

24
 5 Basic graphics functions

>> title ('Subplot 2: Cosine ')

It is possible to let an Axes occupy several grids by letting 𝑝 be a vector.

>> subplot (2 ,2 ,[3 ,4]); % Select the 3rd and 4th grid elements
>> plot(x, y1 , x, y2);
>> title ('Subplot 3 and 4: Sine & Cosine ')

5.6 Other graphics commands

Finally, there are other ways to represent data, and while they will not be elaborated on, it
is worth mentioning. Try running the following script
%% Illustration of a few other graphics functions
% Create data vectors
x = linspace (0 ,2*pi ,30);
y1 = sin(x);

% Open new Figure

figure

% Subplot 1: Discrete sequence plot

subplot (2 ,2 ,1); stem(x, y1);
title('Discrete sequence plot ')

% Subplot 2: Stairstep graph

subplot (2 ,2 ,2); stairs (x, y1);
title('Stairstep graph ')

% Subplot 3: Bar graph

subplot (2 ,2 ,3); bar(x, y1);
title('Bar graph ')

% Subplot 4: Scatter plot

subplot (2 ,2 ,4); scatter (x, y1);
title('Scatter plot ')

25
6 Advanced graphics programming

This chapter will focus on the hierarchical make-up of graphical objects, handles, and
manual construction of Figures. This chapter also introduces the concept of object proper-
ties and ways to configure them. This notion applies to other complex data structures in
Matlab, like the Dynamic System Models, which will be introduced in Chapter 7.

6.1 Graphics objects

Graphics objects are the visual components used by Matlab to create visualizations of
data, where each object plays a specific role in the graphics display.
The relationship of graphics objects is hierarchical and is defined by each object’s parent,
which is in turn a graphics object. Figure 6.1 shows the hierarchical nature of the graphics
objects that will be elaborated on in this text.

Root

Figure

Illustration
Axes
Objects

Chart Primitive
Objects Objects

Figure 6.1: The parent-child relationship between some graphics objects in Matlab1

You can control the behavior and appearance of a particular graphics object by setting its
properties, which will be shown how to do later. First, a brief walkthrough of the important
graphics objects. This list only contains a subset of the most important objects. You can see
the complete list on the Mathworks homepage.

1 The complete image can be found at the Mathworks homepage.

26
 6 Advanced graphics programming

6.1.1 Top-level objects

The Root object is the root of the graphics object tree. Root properties contain information
about the graphics environment and the current state of the graphics system, such as the
monitor configuration, screen size, resolution, etc.
The Figure object is the container for all graphics and represents the display window within
which everything is visualized. Figure properties contain information about the window’s
title, appearance, size, position, etc.
Finally, each Figure should contain at least one Axes object, which serves as the backdrop
upon which lines or other representations of data can be drawn. As already shown in
the previous chapter, the Axes properties control the tick marks, labels, grid lines, and the
relative placement inside the Figure.

6.1.2 Child objects

At the bottom of the tree lie the classes that have no further children.
Chart Objects are created by higher-level graphics commands, like those shown in the
previous chapter. Examples of Chart objects are Line (created with plot), Scatter, Stair,
Stem, etc.
Primitive Objects are simpler objects which make up the above-mentioned Chart objects,
such as Primitive Line objects, Polygons, Images, Text, etc.
Lastly, Illustration Objects are sometimes used to give the reader help interpreting the
graphics. The Legend is an example of such an object.

6.2 Object handles

By default, when graphics objects are created, they will not be directly accessible through
the Workspace, but you can choose to return a handle upon creation. Try the following
commands
>> clear all , close all , clc % Clean slate
>> x = 1:10; % Declare vector
>> p = plot(x, x.^2); % Return chart handle 'p'

The variable p now contains a way to access the Chart Line object created by plot. Using
this handle, you can now directly access the properties of the chart.
The handle only points to the Line object, but as shown in the tree in Figure 6.1, the actual
data representation consists of several other objects. As we have already seen, when calling
plot directly, it will implicitly set up a Figure and an Axes in which to place the Line. It
is still possible to access these directly, but to do so we first need to produce a handle for
them as well.

27
 6 Advanced graphics programming

6.2.1 Get current

One such way is to use the get-current-axes and get-current-figure functions, which will
return handles to the last used Axes and Figure objects, respectively.
>> a = gca; % Get handle to current axes
>> f = gcf; % Get handle to current figure

Check the Workspace and you will see that you now have three graphics objects available
as variables.

6.3 Object properties

Now let’s look at some of the things you can do by accessing the properties directly. Type
>> p

to get a list of some of the properties available. You should see something like this
p =

Line with properties :

Color : [0 0.4470 0.7410]

LineStyle : '-'
LineWidth : 0.5000
Marker : 'none '
MarkerSize : 6
MarkerFaceColor : 'none '
XData : [1 2 3 4 5 6 7 8 9 10]
YData : [1 4 9 16 25 36 49 64 81 100]
ZData : [1 x0 double ]

6.3.1 Setters and getters

To interact with object properties, you can use the default setter and getter methods, where
you specify the property name as a string.
>> X = get(p,'XData ') % Return the x-data as a vector
>> set(p,'Color ','red ') % Set the line colour to red

6.3.2 Dot notation

Since Matlab R2014b, you can use dot notation to access properties. You can revisit
Section 3.6 for more examples of dot notation.
>> Y = p. YData % Return the y-data as a vector
>> p. LineWidth = 2.0 % Increase the line width

28
 6 Advanced graphics programming

6.3.3 Defining properties at declaration

To make objects with custom properties, it is possible to define them at the same time that
the object is created. To define additional properties, you need to specify the key-value
pair, similar to using the setter method shown above.
>> p = plot(x, x.^2 , 'Color ','red ','LineWidth ',1,'LineStyle ','--');

6.4 Creating a publication-quality Figure

Finally, let’s bring this chapter to a close, by studying an example of how to produce quality
figures in Matlab, that integrate smoothly in a document typeset with LATEX. This section
will give an example of what you can do when explicitly controlling the graphics functions.
Figure 6.2 shows the Figure we will create, and the code will be explained, step by step.
The entire code can be found in Appendix 9.8.

Figure 6.2: Vector image created using LATEX interpreting for text fields

In particular, this will cover how to

i) use custom colors, line styles, and thickness for data,
ii) take full control of the figure’s size and the internal axes placement,
iii) customize the axes ticks, labels, and grid lines,
iv) use the LATEX interpreter for typesetting equations,
v) print figure to file, using a vector format.
It is the author’s strongly held view that anyone with the ambitions of producing high-
quality technical reports, research papers, or articles, use this method of generating Figures
if working with data in Matlab.

29
 6 Advanced graphics programming

Preamble
The first thing is to create some data to plot. As the Figure suggests, we are looking at
an exponential smoothing curve and its derivative. We also take the time to define some
custom colors, these are borrowed from the Python package matplotlib.
clear all , close all , clc % Clean slate

% Custom colors
c1 = [31 , 119 , 180]/255; % #1 f77b4
c2 = [255 , 127 , 14]/255; % # ff7f0e

% Data to plot
T = 1; % Time constant (sec)
t0 = 0; tf = 6*T;
t = linspace (t0 ,tf ,50); % Time vector for simulation
f1 = 1 - exp(-t./T); % Exponential smoothing
f2 = exp(-t./T)/T; % Exponential decay
tT = [t0 ,t0+T]; yT = [0 ,1]; % Help - lines for gradient at t=0

Figure dimensions and placement

To make sure the font sizes remain consistent throughout, we will define the Figure in
terms of absolute size. LATEX knows multiple ways to measure lengths, and you can find
out exactly how much vertical space a page fits by using \the\textwidth. See this question
on Stack Exchange for an example how to transform between them.
This page has a \textwidth of 15.27274 cm, so I will let the Figure be 12 cm. The golden
ratio will provide an aesthetic aspect ratio. It is also helpful to place the figure in the center
of the screen.
% Figure size in cm
width = 12;
height = width * 2/(1+ sqrt (5));
pos = [0, 0, width , height ];

% Screen placement
set(groot ,'Units ','centimeters ')
scrn = get(groot ,'ScreenSize ') - pos;
pos (1:2) = scrn (3:4) /2;

Create Figure programmatically

Then we create our Figure, using the dimensions and placement calculated earlier. To use
centimeters, you have to change the units used for the Figure object. All other Objects are
added incrementally, and their Parent explicitly specified, so as to not create undefined
behavior.

30
 6 Advanced graphics programming

Observe that the Text objects, such as the labels, line numbering, and the legend are
instructed to use the LATEX interpreter. This is what makes it possible to typeset equations
in the text fields.
% Create Figure
fig = figure ('Units ','centimeters ','Position ',pos);

% Add Axes
ax = axes('Parent ',fig ,...
'XGrid ','on ','YGrid ','on ','XColor ','k','YColor ','k',...
'XMinorTick ','on ','YMinorTick ','on ','YLim ' ,[0 ,1.2] ,...
'TickLabelInterpreter ','latex ','FontSize ' ,10);

% Set axis labels

set([Link] ,'String ','Time (sec)','Interpreter ','latex ');
set([Link] ,'String ','Amplitude ','Interpreter ','latex ');

% Add Lines
ln (1) = line(ax ,t,f1 ,'Color ',c1 ,'LineWidth ' ,1.0,'LineStyle ','-');
ln (2) = line(ax ,t,f2 ,'Color ',c2 ,'LineWidth ' ,1.0,'LineStyle ','-');
ln (3) = line(ax ,tT ,yT ,'Color ','k','LineStyle ','--');

% Add Legend
leg = legend ('Parent ',fig ,...
'String ',{'$f(x) = 1-e^{-t/T}$','$f''(x) = 1/T\,e^{-t/T}$'},...
'Interpreter ','latex ','FontSize ' ,10,'Location ','southeast ');

Print Figure to file

Finally, we fit the Axes inside the Figure, which we export as a vector image.
% Resize Axes to fill Figure
p = get(ax ,'TightInset '); % Read only!
set(ax ,'Position ',[p(1) ,p(2) ,1-p(1) -p(3) ,1-p(2) -p(4) ]);

% Print Figure to file , using the vector format

print(fig ,'-depsc ','-painters ','PublicationQualityFigure .eps ');

31
7 The Control System Toolbox

The Matlab IDE is designed with high extensibility, and additional features come bundled
in “toolboxes” that can be installed alongside the core system. Up until now, we have been
only using functionality contained in the core toolbox. This chapter focuses on some of the
functions, data structures, and apps contained in the Control System Toolbox. In particular,
only continuous-time, single-input-single-output systems will be addressed in this chapter.
For more advanced topics, see Chapter 8.
Before you continue, make sure you have the Control System Toolbox installed. You can
list all installed toolboxes with the following command
>> ver

Now would be a good time to also check for the other packages listed in Chapter1.

7.1 Dynamic system models

The Control System Toolbox introduces new data structures for representing linear time-
invariant (LTI) systems, some of which are shown in Table 7.1.

Table 7.1: Some of the LTI models in the Control System Toolbox

tf Transfer function model

zpk Zero-pole-gain model
ss State-space model
pid PID controller in parallel form

These numeric LTI models belong to the class of objects known as dynamic system models.
Let’s look at how to create these datatypes, some of their properties, and how to transform
between them. Please note that this chapter will only cover parametric, continuous-time
system models.
Before we continue, let’s wipe the slate clean.
>> clear all , close all , clc

32
 7 The Control System Toolbox

7.1.1 Transfer function models

One of the most common ways to represent LTI systems within systems and control theory
is the transfer function representation. It consists of a fraction of two polynomials of the
complex variable 𝑠, generally written in the form

𝐵(𝑠) 𝑏 𝑚 𝑠 𝑚 + 𝑏 𝑚−1 𝑠 𝑚−1 + · · · + 𝑏1 𝑠 + 𝑏0

𝐺(𝑠) = = .
𝐴(𝑠) 𝑎 𝑛 𝑠 𝑛 + 𝑎 𝑛−1 𝑠 𝑛−1 + · · · + 𝑎 1 𝑠 + 𝑎 0

Transfer functions are used to describe single-input-single-output (SISO) systems, and

mathematically, they describe the dynamics from the input quantity to the output quantity.
For more details on transfer function representations, see your control theory textbook.
You can create a transfer function model in multiple ways. The most straightforward
way involves specifying the coefficients of the numerator and denominator polynomials.
Consider the transfer function
2𝑠 +6
𝐺(𝑠) = 2 .
𝑠 + 8 𝑠 + 15

You can create a model of this using the following commands

>> b = [2, 6]; % Numerator polynomial
>> a = [1, 8, 15]; % Denominator polynomial
>> tf1 = tf(b, a) % Transfer function

Alternatively, you can define the Laplace variable 𝑠 using a reserved keyword, and use this
to construct the rational expression in a (sometimes) more intuitive way.
>> s = tf('s'); % Laplace variable
>> tf2 = (2*s + 6) /(s^2 + 8*s + 15) % Transfer function

7.1.2 Zero-pole-gain models

A crucial aspect of control theory lies in the analysis of the roots of the transfer function
polynomials, which will be elaborated on in Section 7.3.3. The roots of the numerator
polynomial are called the zeros of the system, and the roots of the denominator polynomial
are its poles. Oftentimes, it is convenient to factor the polynomials of the transfer function,
to more easily identify the poles and zeros. Any scaling factor that appears after the
factorization is separated as the system’s gain. The alternative representation can be written
as
(𝑠 − 𝑧 1 ) (𝑠 − 𝑧 2 ) · · · (𝑠 − 𝑧 𝑚 )
𝐺(𝑠) = 𝐾 .
(𝑠 − 𝑝 1 ) (𝑠 − 𝑝 2 ) · · · (𝑠 − 𝑝 𝑛 )

For this purpose, the zpk model can be used. Such a model can be created by specifying
the appropriate arguments, but it is also possible to transform between representations by
using a system model as the input instead. This is the simplest way to transform between
representations and works for all dynamic system models (if such a transform exists). This
method of converting between datatypes will be referred to in this text as casting.

33
 7 The Control System Toolbox

>> zpk1 = zpk(tf1) % Cast TF to ZPK model

The same system can be constructed by specifying the poles and zeros explicitly as vectors.
>> z = -3; % Zero at s = -3
>> p = [-3, -5]; % Poles at s = -3, -5
>> k = 2; % Gain of 2
>> zpk2 = zpk(z,p,k) % Explicitly create ZPK model

With real numbers, a term that appears in both the numerator and denominator in the
fraction can be simplified away. This is also true when considering transfer functions,
and the resulting simplified model is the minimal realization of the original system. In our
example, minreal will perform the pole-zero pair cancellation at 𝑠 = −3.
>> zpk3 = minreal (zpk2) % Minimal realisation

When working with transfer functions, it is often a good habit to always apply a minreal
to the result before studying it further.

7.1.3 State-space models

Another way to represent dynamic systems is the linear state-space representation. This is
a generalized, first-order ordinary differential equation (ODE) of the form

𝑥(𝑡)
¤ = 𝐴 𝑥(𝑡) + 𝐵 𝑢(𝑡)
𝑦(𝑡) = 𝐶 𝑥(𝑡) + 𝐷 𝑢(𝑡),

where 𝑡 is the independent variable (time), and the dependent, vector-valued functions 𝑥(𝑡)
and 𝑢(𝑡) are denoted the state and input vectors, respectively. The vector field 𝑥(𝑡)
¤ defines
the gradient of the state vector at each point in the controlled state-space. In this particular
class of ODEs, the vector field is linear and time-invariant, and so can be represented by
the multiplications with the system (or state) matrix 𝐴, and control (or input) matrix 𝐵.
The function 𝑦(𝑡) defines the output of the system, as some linear function of the state
and input vectors. This is written as the multiplication with the output matrix 𝐶, and
the feedthrough matrix 𝐷 (seldom present). This chapter covers only SISO systems, so the
input and output matrices are vectors. For details on the general MIMO case, see Chapter 8.
A state-space model can be created in Matlab with the ss command. We can cast any of
the models we have seen to a state-space representation, just like earlier. You can specify
additional keywords as strings to elicit desired behavior—see the documentation. Try the
following examples
>> ss1 = ss(tf1) % Cast TF to SS model
>> ss2 = ss(tf1 ,'minimal ') % Minimal realisation

34
 7 The Control System Toolbox

Let’s consider another system, shown below. For the sake of simplifying the notation, the
dependent variable 𝑡 is dropped from future examples.
   
0 1 0
𝑥¤ = 𝑥+ 𝑢
−1 −1 1
𝑦= 1 0 𝑥
 

You can create a state-space model of this by explicitly defining the numerical matrices,
and passing them in as individual arguments.
>> A = [0 ,1; -1 , -1]; % System matrix
>> B = [0;1]; % Control matrix
>> C = [1 ,0]; % Output matrix
>> D = 0; % Feedthrough matrix
>> ss3 = ss(A,B,C,D) % Create SS model

7.1.4 PID controller models

The most classic control architecture within the field of feedback control is the PID con-
troller. It gets its name from the fact that it applies a proportional (P) control signal, as well
as integral (I) and derivative (D) action, the latter usually run through a low-pass filter.
There are commonly two ways to represent a PID controller as a transfer function. The
most common is called the parallel form, and considers the gains of each term individually,
{𝐾 p , 𝐾i , 𝐾d }. The second is called the standard form and uses the same gain 𝐾 p for all
terms, but individual time constants for the integral and derivative terms, {𝑇i , 𝑇d }. If
present, the low-pass filter on the derivative is specified either as the time constant 𝑇f , or
the dimensionless filter divisor 𝑁.
The respective transfer functions are written as
!
𝐾i 𝐾d 𝑠 1 𝑇d 𝑠
𝐶(𝑠) = 𝐾 p + + 𝐶(𝑠) = 𝐾 p 1+ +
𝑠 1 + 𝑇f 𝑠 𝑇i 𝑠 1 + 𝑇d 𝑠
𝑁

where the control parameters are related through 𝐾 i = 𝐾 p /𝑇i , 𝐾d = 𝐾 p 𝑇d and 𝑇f = 𝑇d /𝑁.
While the parallel form is more general and easier to grasp, it should be pointed out
that in the standard form, each control parameter has a clear physical interpretation and
well-defined unit of measurement.
PID controller models can be created with the commands pid and pidstd, depending
on what representation you would like. It is possible to cast existing dynamic system
models to PID form, if such a representation is possible. Any PID model can be cast to the
representations presented earlier.
When creating new models, the trailing arguments can be omitted to create simple P or PI
controllers.

35
 7 The Control System Toolbox

>> Kp = 1; % Proportional gain

>> Ki = 2; % Integral gain
>> pid1 = pid(Kp , Ki) % PI controller with no derivative term

They won’t be elaborated on further in this text, but 2 degrees-of-freedom (2-DOF) PID
controllers are also present in the Control System Toolbox, available with the commands
pid2 and pidstd2.

7.2 Interconnecting dynamic systems

Now that we have seen how to create the individual system representations, it is time to
look at how to connect them, and how to get the resulting numeric representation of the
interconnected systems. First, let’s create some new systems
>> sys1 = zpk ([] , -1 ,1); % ZPK model with unity gain , pole at -1
>> sys2 = zpk ([] , -3 ,3); % ZPK model with unity gain , pole at -3

Figure 7.1 shows three possible ways to connect systems. There are commands in Matlab
for this, but it is also possible to perform arithmetic operations, such as the addition and
multiplication of dynamic system models.

sys1 Σ sys1
−
Σ sys1 sys2

sys2 sys2

(a) Parallel connection (b) Series connection (c) Feedback connection

Figure 7.1: Three general modes of interconnecting system models

7.2.1 Parallel connection

In a parallel connection, the resulting output is simply the sum of the individual system
outputs, under the same input.
>> sys_par1 = parallel (sys1 , sys2) % Systems connected in parallel
>> sys_par2 = sys1 + sys2 % Same as addition

7.2.2 Series connection

In a series connection, the output of the first system is used as the input to the second system.
>> sys_ser1 = series (sys1 , sys2) % Systems connected in series
>> sys_ser2 = sys2 * sys1 % Same as multiplication

Note: The order of the multiplication is important in the multi-variable cases since matrix
multiplication is not commutative in general! The correct order is determined by the reverse
order in which the blocks appear in the diagram.

36
 7 The Control System Toolbox

7.2.3 Feedback connection

When creating a feedback connection, the output of the first system is used as the input to the
second system, upon which the output is connected back to the input of the first system. The
default is to use negative feedback since positive feedback generally results in an unstable
interconnected system.
>> sys_fb1 = feedback (sys1 , sys2) % Systems in a feedback loop
>> sys_fb2 = feedback (sys1 , sys2 , -1) % Explicit negative feedback

It should be mentioned that while not as simple, it is possible to arrive at the same result
with arithmetic operations.
>> one = eye(size(sys1*sys2)); % Identity matrix
>> sys_fb3 = (one+sys1*sys2) \ sys1 % Equivalent equation

To generalize the arithmetic expression for multi-variable systems, the identity matrix is
used instead of just the scalar 1. Also, note the backslash division!

7.3 Linear analysis

Now that we have familiarised ourselves with working with system models, let’s look at
how to analyze the dynamics of the linear systems they represent. This section presents
some tools from the Control System Toolbox that belong to the area of linear analysis, the
ones of most interest being those that produce graphical representations of the system
characteristics. A list of some of the common graphical analysis commands can be seen in
Table 7.2, and a few examples of the resulting Figures can be seen in Figure 7.2.

Table 7.2: Plot commands for linear analysis of dynamic system models

step Step response plot Time response

stepinfo Step-response characteristics Time response
impulse Impulse response plot Time response
lsim Arbitrary input response plot Time response
bode Bode diagram Frequency response
nyquist Nyquist diagram Frequency response
dcgain Steady-state system gain Frequency response
pzmap Pole-zero-map Pole-zero analysis
damp Natural frequency and damping Pole-zero analysis

We can divide the analysis into three domains of study, but it should be stressed that they
are simply different viewpoints of the same fundamental question; “How can we expect
the modeled system to respond to inputs?”
For the sake of the coming analysis, let’s consider the second-order transfer function

37
 7 The Control System Toolbox

Step Response
Bode Diagram
1.6
10

1.4 5

0
1.2
-5
1

Magnitude (dB)
Amplitude

-10
0.8
-15

0.6 -20

-25
0.4
-30
0.2
-35
0
0 5 10 15 20 25 30 -40
10-1 100 101
Time (seconds)
Frequency (rad/s)

(a) Step response (b) Bode magnitude diagram

Nyquist Diagram
Pole-Zero Map
3
1.5
0 dB 0.72 0.58 0.44 0.32 0.22 0.1
-2 dB
2 2 dB
1 0.86
Imaginary Axis (seconds-1 )
4 dB -4 dB
1
6 dB -6 dB 0.5 0.96
Imaginary Axis

10 dB -10 dB
20 dB -20 dB
0 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2
0

-1
-0.5 0.96

-2
-1 0.86

-3 0.72 0.58 0.44 0.32 0.22 0.1

-1.5 -1 -0.5 0 0.5 1 1.5 2 -1.5
-1.8 -1.6 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0
Real Axis
Real Axis (seconds -1 )

(c) Nyquist diagram (d) Pole-zero map

Figure 7.2: Example of some of the graphical linear analysis tools in Matlab

>> sys = tf ([1] , [1, 0.4 , 1]);

which is the system used in Figure 7.2.

7.3.1 Time-domain analysis

The most intuitive system analysis is simply to study how the system output develops
over time, given some known input signal. The most common input for the time-domain
analysis is Heaviside’s unit step function 𝜃(𝑡), defined as
(
0, 𝑡 < 0,
𝜃(𝑡) =
1, 𝑡 ≥ 0.

The resulting system response is called the step response, see Figure 7.2a. You can simulate
and plot the step response for any dynamic system model with the step command. If you

38
 7 The Control System Toolbox

ignore the output, it will plot the response in the most recent Axes, otherwise, it will just
save the vector containing the time and outputs, to the Workspace.
>> step(sys) % Step response plot
>> [y,t] = step(sys) % Returns the output and time vectors

For stable LTI systems, the overall behavior is a fast initial rise, followed by convergence
towards some final steady-state value. In the case of our oscillating systems, there is also
some overshoot present. You can extract system characteristics such as rise time, overshoot
in percent, and settling time, by using
>> stepinfo (sys) % Computes step response characteristics

Another common system response is the impulse response, where the input is a unit impulse
at 𝑡 = 0. This can be thought of as an instantaneous transfer of energy into the dynamic
system. Just like the step response, you can choose to produce the plot or to save the output
as vectors.
>> impulse (sys) % Impulse response plot

Additionally, you can test any input signal of your choice by creating the necessary time
and input vectors, and exciting the linear system using lsim.
>> t = 0:0.1:30; % Time from 0 to 30 sec
>> u = asin(sin(t)); % Saw -wave with period 2 PI
>> lsim(sys ,u,t) % Linear simulation plot

7.3.2 Frequency-domain analysis

The other useful domain of analysis for LTI systems is the frequency domain. Recall that
the Laplace variable 𝑠 is a complex number defined as

𝑠 = 𝑗 𝜔 + 𝜎,

where the imaginary part contains some frequency 𝜔, and the real part corresponds to
some factor of decay 𝜎. By letting the real part be zero, we can study how the system
behaves when excited with undamped oscillations at various frequencies. This is called the
frequency response.
The most useful plot for studying the frequency response is the Bode diagram, which shows
the magnitude (dB) and phase shift (degrees) of the system, as a function of the input
frequency (rad/sec). Similar to the previous graphical commands, you can either return
the output vectors to the Workspace or show the plot window.
>> bode(sys) % Bode diagram

Oftentimes, it is sufficient to only show the magnitude plot, like in Figure 7.2b.
>> bodemag (sys) % Only show the magnitude plot

39
 7 The Control System Toolbox

Another way to represent the frequency response is the Nyquist diagram, which plots the
system magnitude and phase shift in the complex number plane, using polar coordinates
and extending the range of frequencies to include negative numbers.
>> nyquist (sys) % Nyquist diagram

An example of the resulting parametric curve can be seen in Figure 7.2c, and is a common
tool for analyzing the robustness of linear control systems.
Lastly, it is possible to extract important numeric values relating to the frequency response
of a system, by using the commands below. The steady-state gain of a system is defined as
the magnitude of the system when excited with a frequency of zero, hence why it’s also
referred to as the DC gain, from the field of electrical engineering.
It is equivalent to setting 𝑠 = 0 in the transfer function and can be found with the command
>> K0 = dcgain (sys) % Steady - state gain

Another important system characteristic is the bandwidth, which is the frequency where the
gain drops below =3 dB of its steady-state value.
>> B = bandwidth (sys) % Bandwidth of the system

Any frequencies above the bandwidth are commonly considered strongly attenuated.

7.3.3 Pole-zero analysis

Without considering any inputs explicitly, it is also possible to draw conclusions regarding
the system behavior by just looking at the roots of the transfer function polynomials. One
such characteristic is stability. Quoting from the MathWorks documentation;
Stability is a standard requirement for control systems to avoid loss of control
and damage to equipment. For linear feedback systems, stability can be assessed
by looking at the poles of the closed-loop transfer function.
A dynamic system modeled as a rational transfer function of order 𝑛 is stable if and only if
all the system poles 𝑝 have a negative real part;

Re (𝑝 𝑖 ) < 0, 𝑖 = 1, 2, . . . , 𝑛.

You can determine the pole and zero locations of a SISO system model using the commands
>> p = pole(sys) % Poles of the system
>> z = zero(sys) % Zeros of the system

A useful way to graphically represent the poles and zeros of a transfer function is by
plotting their locations on the complex 𝑠-plane, whose 𝑥 and 𝑦-axes represent the real and
imaginary parts, respectively. It is common to mark a zero location with a circle (◦), and a
pole location with a cross (×). Such a figure is known as pole-zero map, see Figure 7.2d.
>> pzmap (sys) % Pole -zero map

40
 7 The Control System Toolbox

To ensure system stability, no poles can be to the right of the line drawn by the imaginary
axis. Additionally, other conclusions regarding the dynamics of a stable system model can
be drawn by studying the pole locations more carefully. Since all the coefficients of the
transfer function polynomials are real, the poles and zeros must be either purely real or
appear in complex conjugate pairs.
• A purely real pole corresponds to an asymptote in the time domain, where the system
output converges exponentially to some steady-state value. The rate is determined
by the time constant, which is the inverse of the real pole.
• Complex-conjugated pole pairs represent oscillations in the time-domain behavior.
The magnitude and angle of the complex phasor can tell you about the natural frequency
of the oscillations, as well as the damping ratio.
You can determine the numerical values using the following command
>> damp(sys) % Natural frequency and damping of the system

For more theory regarding conjugated poles and oscillations, consult the literature from
your basic control theory course.

Exercise
Consider the stable transfer function 𝐺(𝑠) in a feedback control-loop with the controller
𝐶(𝑠) shown in Figure 7.3.

𝑟(𝑡) 𝑦(𝑡)
Σ 𝐶(𝑠) 𝐺(𝑠)
−

Figure 7.3: Block diagram of the feedback-controlled system

The transfer functions are given below

100 4 𝑠 2 + 40 𝑠 + 400
𝐺(𝑠) = 𝐶(𝑠) =
𝑠 + 10 𝑠 + 100
2 𝑠 2 + 40 𝑠

Do the following in a new Matlab script, and try to answer the questions below.
1. Create the transfer function models for 𝐺(𝑠) and 𝐶(𝑠).
2. Create the resulting closed-loop system from 𝑟(𝑡) to 𝑦(𝑡), shown in Figure 7.3.
3. Investigate the closed-loop system poles and zeros. Notice anything particular?
Hint: It might be more clear if you represent it as a zero-pole-gain model.
4. Now try to determine the minimal realization of the closed-loop system. What
changed?

41
 7 The Control System Toolbox

5. The controller is actually just a PID controller. What are its parameters?
6. Compare the time-domain properties of the controlled system (such as rise time,
overshoot, etc.) of the uncontrolled system 𝐺(𝑠). Can you see some correlation if you
then compare the frequency-domain properties such as poles, damping, etc.?

42
8 Dynamic systems – Advanced topics

This chapter covers some of the more advanced aspects of modeling dynamic systems in
Matlab. In particular, this chapter serves to extend some of the definitions provided in
the previous chapter, to cover more advanced topics such as multi-variable, and discrete-
time systems. Lastly, it will present alternative commands for working further with the
graphical representations shown in the previous chapter.
To provide a preface, the first section covers how to manipulate dynamic system models,
and access their internal properties. For this reason, it is suggested that you first work
through Chapter 6, if you have not already.
Please note that this chapter involves some minor abuse of notation, where time-domain
signals and their transformed counterpart are used interchangeably. For convenience, the
time dependency and the Laplace variable may also be dropped from the notation.

8.1 Model objects

As alluded to at the beginning of the previous chapter, dynamic system models are objects
and contain properties that can be read or explicitly defined. Let’s create a new transfer
function using these methods. While not the most convenient method, it helps to illustrate
the general procedure of accessing the object properties.
Begin by creating a new, empty transfer function.
>> clear all , close all , clc % Clean slate
>> tf1 = tf % Empty transfer function

The simplest way to access the properties is by using the dot-notation. First, we give our
model a name.
>> [Link] = 'Low -pass filter ';

Alternatively, you we use the set and get commands. We can define the polynomials by
setting the name-value pair explicitly.
>> set(tf1 , 'Numerator ', 10, 'Denominator ', [1, 10]);

The default time unit is seconds. Confirm that with the appropriate get command.
>> get(tf1 ,'TimeUnit ') % Read the time unit

43
 8 Dynamic systems – Advanced topics

The raw numeric representations can also be extracted from LTI models with the commands
tfdata, ssdata, etc. depending on the type of model.
>> [num , den] = tfdata (tf1) % Get the polynomial coefficients

8.2 Multi-variable systems

You will notice that the coefficients are contained inside a 1 × 1 cell array (for details on
cells, see Section 3.5). The reason for this is that transfer functions (like almost everything
else in Matlab) can be indexed, and used to represent multi-variable systems.
In general, multiple-input-multiple-output (MIMO) systems are characterized by having
vector-valued inputs and/or outputs. The SISO transfer function is simply a scalar (1 × 1)
instance of a more general transfer matrix object.
Before we consider transfer matrices, the more intuitive approach to understanding MIMO
systems is by first considering the state-space representation, since it is simple to extend to
the MIMO case.

8.2.1 State-space models

The dimensions of a state-space model denote the lengths of the vectors which make up
the model. For linear systems, the number of internal states, inputs, and outputs are
determined by the state matrix 𝐴, the input matrix 𝐵, and the output matrix 𝐶, respectively.
The matrix dimensions must agree since the system is solved numerically using matrix
multiplication. Given 𝑝 states, 𝑞 inputs, and 𝑟 outputs, the requirements on vector and
matrix dimensions are listed in Table 8.1.

Table 8.1: State-space dimension requirements

𝑥 : 𝑝×1 𝐴:𝑝×𝑝
𝑥¤ : 𝑝 × 1 𝐵:𝑝×𝑞
𝑢 : 𝑞×1 𝐶 :𝑟×𝑝
𝑦 : 𝑟×1 𝐷:𝑟×𝑞

Physical systems rarely have feedthrough characteristics, so very often the 𝐷-matrix is just
a zero matrix. In newer versions of Matlab you can simply specify 0 and the size will be
determined automatically. For older versions, the following dynamic sizing can be used
D = zeros (size(C ,1) , size(B ,2));

Let’s re-consider the system model from Chapter 7, but as a multi-variable system.
>> A = [0 ,1; -1 , -1]; % System matrix
>> B = [0;1]; % Control matrix
>> C = [1 ,0;0 ,1]; % Output matrix
>> ss1 = ss(A,B,C ,0) % Create SS model

44
 8 Dynamic systems – Advanced topics

Setting the identity matrix as the output matrix ensures that the entire state vector is
returned as the output (𝑦 = 𝑥) since a vector multiplied by the identity matrix is unchanged.
By specifying the index in Matlab, you can specify which output you want, try
>> ss1 (1) % Model of the dynamics from u to y1
>> pole(ss1 (2)) % Poles of the model from u to y2
>> step(ss1) % Step response from both outputs

Next, let’s make the system multi-input as well. We can change the control matrix 𝐵 to be
multi-variable without re-initializing the model, but this requires setting the 𝐷 matrix too,
to avoid matrix inconsistencies
>> set(ss1 , 'B', eye (2) , 'D', zeros (2)) % Vector - valued input

Setting the control matrix 𝐵 to the identity matrix means that the system now takes the
vector-valued input 𝑢 = [𝑢1 , 𝑢2 ]T , where 𝑢1 acts directly on the state 𝑥1 , and 𝑢2 on 𝑥 2 .
When working with a MIMO system, you can choose to study individual input-output
dynamics by providing two indices, the first for the output, and the second for the input.
This ordering is no coincidence, as will be shown in the next section. The syntax for
indexing is just like with regular matrices, see Section 3.3.4.
>> ss1 (1 ,1) % Model of the dynamics from u1 to y1
>> step(ss1 (: ,2)) % Unit step at u2 , show both outputs
>> step(ss1) % Unit step at each input , respectively

By providing the signals for both inputs explicitly, it is possible to simulate inputs from
both inputs simultaneously
>> t = linspace (0 ,10 ,10000) % Time vector
>> u = ones (2 ,10000); % Unit step at both inputs
>> lsim(ss1 ,u,t) % Linear simulation

8.2.2 Transfer matrices

Another way of modeling MIMO systems is using transfer matrices, where each element of
the matrix contains a single SISO transfer function. Consider some dynamic system with
𝑞 ∈ ℕ inputs and 𝑟 ∈ ℕ outputs. The resulting system of equations can be written as

𝑦 = 𝐺 𝑢,

where 𝑦 is the output vector, and 𝑢 is the input vector. The transfer matrix 𝐺 assumes the
dimensions 𝑟 × 𝑞, and the form

𝑔1,1 (𝑠) 𝑔1,2 (𝑠) · · · 𝑔1,𝑞 (𝑠)

𝑔2,1 (𝑠) 𝑔2,2 (𝑠) · · · 𝑔2,𝑞 (𝑠)
 
𝐺(𝑠) =  . .. .. ..  .
.

 . . . . 
𝑔 (𝑠) 𝑔 (𝑠) · · · 𝑔 (𝑠)
 𝑟,1 𝑟,2 𝑟,𝑞 

45
 8 Dynamic systems – Advanced topics

It becomes clear that the transfer function from input 𝑢 𝑗 to output 𝑦 𝑖 is given by 𝑔 𝑖,𝑗 , which
denotes the entry in the 𝑖 th row and the 𝑗 th column of the transfer matrix. For this reason,
it becomes natural to index all MIMO systems in Matlab this way, which is why indexing
for state-space systems follows this convention too.
You can cast a state-space model to a transfer matrix just like with scalar systems
>> tf1 = tf(ss1) % Cast to transfer matrix
>> zp1 = zpk(ss1) % Also works in ZPK form

You will see that you can do the same things with the transfer matrix, as with the state-space
model.
>> tf1 (1 ,1) % Transfer function from u1 to y1
>> step(tf1 (: ,2)) % Unit step at u2 , show both outputs
>> step(zp1) % Unit step at each input , respectively

8.3 Discrete-time systems

So far, all models mentioned have been in continuous-time, where signals can assume a
particular value for an infinitesimally short amount of time. It should be mentioned that
most mechatronic applications are computer-controlled systems, where continuous signals
are converted to digital representations, resulting in discrete-time systems.
The fundamental property that separates a discrete system model from its continuous
counterpart is that a discrete system model has a fixed sample time. It is assumed that
the system being modeled only updates its measurements at regular intervals determined
by the sample time, with values assumed constant in between updates. This is called
zero-order hold (ZOH).

8.3.1 The discrete-time domain

In discrete time, continuous differential equations have been replaced with difference equa-
tions, which are used to express a signal as a sequence of samples, with their order de-
termined by the monotonically increasing index 𝑘 ∈ ℤ. Just like the continuous-time
state-space model is a generalization of the ordinary differential equation, the discrete-time
state-space model is a generalization of the difference equation, written as
𝑥[𝑘 + 1] = 𝐴 𝑥[𝑘] + 𝐵 𝑢[𝑘],
𝑦[𝑘] = 𝐶 𝑥[𝑘] + 𝐷 𝑢[𝑘].

The discrete state-space equation describes the value of the state vector at sample 𝑘 + 1 as a
linear function of the state and input vectors of the previous sample 𝑘, defined by its four
system matrices, which take on a form very familiar to that we have already seen.
Similarly, the discrete-time transfer function also closely resembles its continuous counter-
part; a fraction of two polynomials of a complex variable, written as
𝐵(𝑧) 𝑏 𝑚 𝑧 𝑚 + 𝑏 𝑚−1 𝑧 𝑚−1 + · · · + 𝑏1 𝑧 + 𝑏0
𝐺(𝑧) = = .
𝐴(𝑧) 𝑎 𝑛 𝑧 𝑛 + 𝑎 𝑛−1 𝑧 𝑛−1 + · · · + 𝑎1 𝑧 + 𝑎 0

46
 8 Dynamic systems – Advanced topics

Analogously, a discrete transfer matrix describes the transfer from the sequence of inputs
to the sequence of outputs. Note that the complex variable is 𝑧 instead of 𝑠, which serves
as a reminder that discrete-time system models are functions with a different domain.
There are several ways of transforming between the domains, most of which involve some
degree of approximation, where models may lose some of their original expressiveness.
This text is not going to cover any in-depth theory on discretization or 𝑧-transform, for that
purpose, see a textbook in digital control theory1 .

8.3.2 Declaring discrete system models

Discrete-time system models can be declared explicitly, just like in the previous Chapter,
with the sample time provided as an additional argument
>> Kp = 1; % Proportional gain
>> Ki = 2; % Integral gain
>> Ts = 0.01; % Sample time
>> pidd = pid(Kp , Ki , 0, Ts) % Discrete PI controller

It is important to note that in the discrete-time, different rules apply. For instance, a pole in
the right-half plane will produce a stable system, as long as its magnitude is less than one.
>> z = 0; % Zeros
>> p = 0.5; % Poles
>> k = 0.5; % Gain
>> Ts = 0.01; % Sample time
>> zpkd = zpk(z, p, k, Ts) % Discrete transfer function

8.3.3 Converting discrete system models

Most of the time, discrete system models are arrived at by discretizing an already existing,
continuous system model. This is done with the c2d command.
>> tfc = tf(1, [1 ,1]); % Transfer function , pole at s = -1
>> Ts = 0.1; % Sample time
>> tfd1 = c2d(tfc , Ts) % Discretise using ZOH method
>> tfd2 = c2d(tfc , Ts , 'zoh ') % Can also be stated explicitly

The default method for discretization method provides an exact discretization in the time
domain, given ZOH acting on the inputs. Another common method is Tustin’s approximation,
which provides a good match in the frequency domain, and allows for frequency pre-
warping to preserve important dynamics at certain frequencies.
>> tfd3 = c2d(tfc , Ts , 'tustin ') % Tustin 's approximation

Finally, it should be mentioned that it is also possible to apply the reverse transformations
to the methods mentioned above, using the c2d command.
>> tfc1 = d2c(tfd1)

1 For a complete overview, see Computer-Controlled Systems: Theory and Design by Åström and Wittenmark.

47
 8 Dynamic systems – Advanced topics

Exercise
Consider the continuous state-space model
   
0 1 0
𝑥(𝑡)
¤ = 𝑥(𝑡) + 𝑢(𝑡),
−1 −1 1
| {z } |{z}
𝐴 𝐵
𝑦(𝑡) = 1 0 𝑥(𝑡).
 
| {z }
𝐶

Using a simple forward Euler approximation

𝑥(𝑡 + 𝑇s ) ≈ 𝑥(𝑡) + 𝑇s · 𝑥(𝑡),

¤

it follows that an approximation of the state-space model is given by

𝑥[𝑘 + 1] = Φ 𝑥[𝑘] + Γ 𝑢[𝑘], Φ ≈ 𝐼 + 𝐴 𝑇s

𝑦[𝑘] = 𝐶 𝑥[𝑘] Γ ≈ 𝐵 𝑇s

Do the following in a new Matlab script, and try to answer the questions below.
1. Create a continuous state-space model of the original system.
2. Create a discrete state-space model of the approximated system with 𝑇s = 0.5 sec.
3. Discretise the original model using c2d with ZOH and the same 𝑇s .
4. Compare the step responses of the three systems, in the same Axes. How good are
they at capturing the model dynamics?
5. Repeat the above for 𝑇s = 0.1 and 𝑇s = 1, what can you observe?
6. Run the script and answer the same questions using the approximation

1
Φ ≈ 𝐼 + 𝐴 𝑇s + (𝐴 𝑇s )2 Γ = 𝐴−1 (Φ − 𝐼) 𝐵
2

7. Optional: You might recognize the approximation as the first two terms in the Maclau-
rin series expansion of a certain function 𝑓 (𝐴 𝑇s ). Use that function instead to deter-
mine Φ (and thus Γ). How does the approximation compare to ZOH now?
Hint: For performing this function on matrices, check
>> help matfun

48
 8 Dynamic systems – Advanced topics

8.4 Customisable plots

Lastly, this section is meant to serve as a complement to Chapter 6, and provide some com-
mands which help facilitate advanced graphics programming for linear analysis figures.
As shown in Section 7.3, the linear analysis functions step, bode, pzmap etc. will not
produce a Figure if any return arguments are specified. Omitting the arguments works
perfectly fine if you are just interested in a quick visualization while developing your code,
however, for the sake of producing quality Figures, another method is preferred.
The Control System Toolbox provides alternative versions of the previously mentioned
commands, which return the handle to the plot. A few are listed in Table 8.2. As you will
notice they all have the added ending “plot”.

Table 8.2: Some plot commands with additional customization options

stepplot Step response plot

lsimplot Arbitrary input response plot
bodeplot Bode diagram
pzplot Pole-zero-map

Once you have the handle, it becomes possible to exercise a high level of control over the
plot’s appearance. It should be mentioned that these plots are specializations of the more
graphical objects shown in Chapter 6, so in addition to modifying the underlying objects
such as the Axes, you may also provide options specific for that particular plot.
For instance, Figure 7.2b can be produced with the following code snippet
sys = tf(1, [1 ,0.4 ,1]); % Transfer function
handle = bodeplot (sys); % Get plot handle
setoptions (handle , 'Grid ','on ', 'PhaseVisible ','off '); % Options

In this case, disabling the PhaseVisible option is used to show only the magnitude, not the
phase shift. For the available arguments, you should consult the MathWorks documentation
for the command in question.
For the actual code used to produce the plots in Figure 7.2, see Appendix 9.8.

49
9 Simulink

Simulink provides a graphical user interface (GUI) for building dynamic system models
as block diagrams, using click-and-drag mouse operations. With this interface, you can
draw the models just as you would with pencil and paper (or as most textbooks depict
them). This chapter explains some of the basic features of the Simulink software package,
examples on how to build a model from scratch, and the integration with the Matlab
engine, including some of the toolboxes and techniques seen in previous chapters.

9.1 Motivation
Simulink encourages you to try things out. You can easily take an existing model and add
to it, or build models from scratch, as will be shown in Section 9.5. Using scopes and
other display blocks, you can see the simulation results while the simulation is running. In
addition, you can change many parameters and see what happens for “what if” exploration.
After you define a model, you can simulate it using a choice of integration methods, either
from the Simulink menus or by entering commands in the Matlab Command Window.
With Simulink, it becomes easier to move beyond idealized linear models, to explore more
realistic nonlinear models, factoring in friction, air resistance, gear slippage, hard stops, and
the other things that describe real-world phenomena. The Simulink modeling environment
supports nonlinear and discontinuous dynamics, modeled in continuous time, discrete
time, or a hybrid of the two. Systems can even be multi-rate, i.e. have different parts that
are sampled or updated at different rates.

9.2 Fundamentals
All models in Simulink are made up of blocks, which may have any number of inputs and/or
outputs. A block can represent a physical component, a small system, or a function.
Additionally, Simulink handles data in three categories, listed below.
• Signals – Block inputs and outputs, computed during simulation.
• States – Internal values, representing the dynamics of the block, computed during
simulation.
• Parameters – Values that affect the behavior of a block, controlled by the user.

50
 9 Simulink

In short, signals represent the connection between blocks’ output(s) and input(s). It is
possible to create a branch in a signal so that the output from one block serves as the
input to several blocks. Signals can be vectors of arbitrary length, and Simulink supports
numerical and logical values.
Similar to Matlab, signal types and dimensions are inferred automatically unless defined
explicitly. When a model is run for the first time it is compiled into executable code, upon
which all parameters and signals are validated. Any mismatched parameter or signal
properties will cause errors, and prevent the simulation from running.
When working in Simulink, there are a few things to keep in mind to try to maintain a
logical and systematic structure to your models:
• Try to work from left to right, with inputs at the very left and outputs at the very
right.
• Label your signals using some consistent convention that is easy to follow.
• Refrain from using hard-coded values in your blocks, use variables defined in the
Matlab workspace.
• Use Subsystems to group blocks, which increases readability and reduces the com-
plexity of your models.

9.2.1 Getting started

You can open a new instance of Simulink in several ways, either by clicking the “Simulink”
button in the Home Toolstrip of the Matlab desktop or by using the command
>> simulink

On the Simulink start page, choose a new “Blank Model”. You should see a blank canvas,
with a Toolstrip at the top, similar to that of Matlab. You will also see that Matlab is still
open, and closing it will prompt Simulink to also close since Simulink is not a standalone
environment.
The canvas can be populated with any number of blocks to create a model. When working
on the canvas, you can use the mouse to
• move blocks by holding the left mouse button,
• copy blocks by holding the right mouse button,
• open block parameters by double-clicking with the left mouse button,
• connect an input and output by holding the left mouse button,
• create a new branch on a signal by holding the right mouse button,
• label existing signals by double-clicking with the left mouse button,
• zoom in/out in the model using the mouse wheel,
• move the canvas by holding the middle mouse button,

51
 9 Simulink

• bring up the context menu by clicking the right mouse button.

9.2.2 Library Browser

You can access the complete list of pre-defined blocks in Simulink by using the Library
Browser, which is available in the Simulation Toolstrip, see Figure 9.1.

Figure 9.1: Simulink Library Browser

Brief descriptions of the most important categories are listed below.

• Commonly Used Blocks – contains a selection of useful blocks present in other groups.
• Continuous / Discrete – contains blocks for building linear time-invariant systems us-
ing transfer functions or state-space form in continuous or discrete time, respectively.
• Discontinuities – contains blocks for representing discontinuous signal properties,
such as Saturation and hysteresis (Relay).
• Math Operations – contains blocks for elementary arithmetic, such as Sum and Prod-
uct, as well as unary functions such as Sqrt, Abs, and Sign.
• Signal Routing – contains blocks for stacking and splitting signals (Mux, Demux) and
switching between different signals during runtime.
• Sinks – contains blocks with no output of their own, for terminating branches and
visualizing (Scope) or exporting the output data (To Workspace).

52
 9 Simulink

• Sources – contains blocks with no input of their own, which provide signal generation
functionality, such as the Pulse Generator, Step, and Sine Wave blocks.

9.2.3 Running a simple example

Let’s begin by adding your first blocks. Drag-and-drop a Sine Wave from the Sources.
Upon releasing the block, you will be prompted to define the first block parameter. In
this case, you can leave it at 1, then double-click the block to open the complete parameter
configuration. Set the frequency to 2*pi rad/sec, which corresponds to 1 Hz. Next, add a
Scope from the Sinks into the canvas, and connect the blocks by drawing an arrow between
the ends using the mouse. The resulting model should look like Figure 9.2.

Figure 9.2: Simple two-block model in Simulink

In the Simulation Toolstrip, in the “Simulate” section, set the Stop Time to 5, and click the
Run button to simulate the model. The progress is shown in the status bar at the bottom
of the window. In this case, the model is not particularly complex and should compile and
execute within a few seconds. Once the execution finishes (and if no errors occurred), the
status will return to Ready, and you may hear a system notification sound. The simulation
results can be seen in the Scope block by double-clicking it; you should see something like
Figure 9.3.

Figure 9.3: Simulation output visualized in the Scope

53
 9 Simulink

The Sine Wave output is displayed in an oscilloscope-style grid, with time along the 𝑥-
axis. It is more than likely that the result does not exactly reflect your expectations. For
instance, the amplitude of the signal never reaches 1, and the curve appears to have several
discontinuous breaks. This will be addressed in the coming section.

9.3 Solver configuration

A crucial element of using Simulink that often is overlooked is the solver and how it is
configured. This text will not delve deep into the field of numerical analysis for ordinary
differential equations, but only touch on some common issues and things to keep in mind1 .
The single most important parameter for solving equations numerically is the step size;
generally, a smaller step size increases the accuracy of the results but also increases the time
required to simulate the system. In Simulink, there are two classes of solvers that can be
used.

9.3.1 Fixed-step methods

This includes all methods where the step size is constant for the entire simulation. Choosing
a fixed-step solver can produce more accurate results, but requires a good understanding
of the system and can cause errors if the step size is chosen incorrectly. Simulink offers
solvers of various orders of accuracy, signified by the trailing number in the name. For
instance, the Runge-Kutta method ode4 has an order of accuracy of four, meaning the total
accumulated error is on the order of 𝑂(ℎ 4 ), where ℎ is the step size. Decreasing the step
size will decrease the bound on the error as defined by the solver, but also increase the
computational burden.
This family of solvers can be used when you have a simple model, or the required computa-
tion for choosing a small step size for the entire simulation is low. Additionally, if you plan
to generate code from your model and run the code on a real-time system with constant
sampling, you should also use a fixed-step solver to reflect the target hardware.

9.3.2 Variable-step methods

Most of the time, a trade-off between accuracy and performance can be achieved, where
the step size is determined dynamically during runtime. This is facilitated by variable-step
methods. Generally, the step size is chosen based on the smoothness of the solution curve
– becoming smaller for periods with fast transients, and larger when the dynamics are
slow or the system is at rest. The flagship solver for ODEs in Matlab and Simulink is
ode45 (order of accuracy 4 to 5), which performs well with most ODE problems and should
generally be your first choice.
However, there are classes of systems of equations for which it may be suitable with another
solver. One such class is the differential-algebraic equation (DAE), which will not be covered
in this text, but often arises when dealing with discontinuities or nonlinearities. Another
important class of systems is the so-called stiff ODE system.

1 For a good summary of key concepts, see Chapters 5 and 6 of Numerical Analysis (2nd ed.) by Timothy Sauer.

54
 9 Simulink

9.3.3 Stiff systems of equations

The stiffness of an ODE is a crucial factor that directly affects the performance of a solver.
When simulating a system with high stiffness, the step size may need to be decreased to an
unacceptably small size in order for a numerical method to give a reliable solution, even in a
region where the solution curve is smooth. This will commonly cause an error in Simulink.
Stiffness is a subtle, difficult, and important concept in the numerical solution of ordinary
differential equations.2 It is hard to tell exactly when a system can be considered stiff, but
a good indicator is that the simulation is unusually slow, or an error message about the
solver decreasing the sample size beyond a tolerated limit. For stable, higher-order linear
systems, a useful metric is the stiffness ratio, defined as the quotient of the smallest and
largest eigenvalues (poles) of the system. To put it simply: the further apart the system
poles, the higher the stiffness.
For dealing with stiff systems in Simulink, there are a few implicit solvers which per-
form more calculations each step, but typically require much fewer steps to complete the
simulation. Knowing when to use which solver is difficult, and often requires some exper-
imentation. Table 9.1 gives a brief description of the methods and recommendations for
when to use them3 .

Table 9.1: General guidelines on when to use each of the different solvers

Solver Problem type Order of accuracy When to use

ode45 Nonstiff Medium Most of the time. This should be the first
solver you try.
ode23 Nonstiff Low For problems with crude tolerances, or in the
presence of moderate stiffness
ode113 Nonstiff Low to high For problems with stringent tolerances, or
when the ODE is expensive to evaluate.
ode15s Stiff Low to medium When ode45 fails or is inefficient and you sus-
pect that the problem is stiff, or when solving
differential-algebraic equations (DAEs).
ode23s Stiff Low When solving stiff systems with constant
mass matrix and known Jacobian, and ode15s
proves inefficient for crude tolerances.
ode23t Moderately stiff Low If the problem is only moderately stiff and you
need a solution without numerical damping,
or when solving DAEs.
ode23tb Stiff Low When solving stiff systems using crude toler-
ances, and ode15s proves inefficient.
ode15i Fully implicit Low For fully implicit problems 𝑓 (𝑡, 𝑦, 𝑦 ′) = 0 and
differential algebraic equations of index 1.

2 For more details, see the article Stiff Differential Equations by Cleve Moler on the MathWorks homepage.
3 Adopted from the page Choose an ODE Solver on the MathWorks homepage.

55
 9 Simulink

9.3.4 Discrete solvers

Additionally, for systems with no continuous states, there are discrete solvers with both
fixed and variable step sizes. In discrete time, the solution curve does not need to be solved
numerically, since the nature of discrete systems allows the state evolution to be uniquely
determined from the difference equation alone. These solvers can only be used when the
entire model is in discrete time, causing an error otherwise.

9.3.5 Automatic selection by Simulink

Unless configured explicitly, the default behavior of a Simulink model is to try to automat-
ically determine a variable-step solver that best suits the dynamics of the model. This is
done using the heuristics seen in Table 9.2, and gives some insight into where to begin4 .

Table 9.2: Auto solver heuristics

Model has continuous states

has only discrete states
Solver is a system of ODEs contains DAEs

Fixed-step FixedStepDiscrete ode3 ode14x

is nonstiff is stiff
Variable-step VariableStepDiscrete ode23t
ode45 ode15s

This also brings us back to the original issue, the poor resolution of the Sine Wave model.
If you go back to your model, you will likely see that the status bar is displaying
auto(VariableStepDiscrete)
in the bottom right. Knowing what we know now, it makes sense that this is the solver it
chose – the Sine Wave has no continuous states since it is just a Source.
If you click the text, you will also see that it has automatically determined the Max step size
to be 0.1, which explains why the curve has such low resolution. Click the cogwheel icon
(see Figure 9.4), and it will take you to the solver settings. Alternatively, you can get there
from the Modeling Toolstrip, under the Setup tab.

Figure 9.4: Bringing up the solver settings from the status bar

4 Adopted from the page Choose a Fixed-Step Solver on the MathWorks homepage.

56
 9 Simulink

If you expand the Solver details you will get access to all the solver parameters. Try setting
the Max step size to 0.001 and running the simulation again. Hopefully, the curve should
look more like a sinusoidal wave than before.
Hint: You can run the simulation directly from the Scope, by opening it again and using
the Run button in the Scope toolbar.

9.4 Studying a demo model

Now that we have covered the fundamentals of the Simulink environment, let’s look at how
a dynamic system is modeled in Simulink. A good place to start exploring the intricacies of
dynamic systems modeling is to try one of the included demo models. Close the untitled
Sine Wave model, and type the following command in the Matlab Command Window
>> sldemo_househeat

This will start a new instance of Simulink and open the system model shown in Figure 9.5.
blower
Terr cmd
On/Off
HeatFlow 1/s
Fahrenheit HeatCost
Thermostat Troom
to Celsius Cost
Set Point
Calculator PlotResults
Heater
F2C
C2F

House
F2C Tindoors
Tout
Avg Outdoor
Temp

Daily Temp
Variation Toutdoors
Temperatures

Copyright 1990-2012 The MathWorks, Inc.

Figure 9.5: Thermal Model of a House, a demo model from MathWorks

The demo models the (simplified) thermodynamics of a house, and the overall dynamics
can be described as follows.
1. The Thermostat controls the Heater, which is switched on or off based on the difference
between the reference (Set Point) and measured temperature.
2. The temperature of the House is affected by the output from the Heater, as well as
the outside temperature, which fluctuates over the course of a day. This is simulated
as a sine wave superimposed on a daily average temperature.
3. Finally, the temperature of the House is measured and fed back into the Thermostat,
which closes the loop on this simplified feedback control system.

57
 9 Simulink

9.4.1 Coordinate transforms

It should be mentioned that all external temperatures are measured in degrees Fahrenheit,
whereas the internal dynamics are specified in degrees Celsius. This is achieved using the
necessary C2F and F2C blocks, which perform the corresponding transformations between
the two measurement systems.
This kind of “wrapped” architecture is common in other applications too, where the trans-
form may be between other frames of reference or coordinate systems.

9.4.2 Subsystems
The model uses subsystems to simplify the model diagram and create reusable systems. A
subsystem is a grouping of blocks, implemented by the Subsystem block. If you double-
click a Subsystem block it will show the contents of the underlying blocks. The oval blocks
represent the input and output signals, and the number inside the block indicates the port
number of that signal.
This model contains three main subsystems, shown in Figure 9.6; representing the Ther-
mostat (9.6a), Heater (9.6b) and House (9.6c), respectively.

1
Thermostat Command (On/Off)
On/Off

1 1 1
2
HeatFlow
Terr Blower Troom Heater Switch
switch HeatGain
(a) Thermostat Subsystem
Heater Air Temperature

(b) Heater Subsystem

1 1/s 1
Heater Room Temp
QDot Troom
In 1/Mc

Heat Losses
2
Outdoor Temp
Tout

(c) Thermodynamic model for the house

Figure 9.6: Subsystems making up the demo Thermal Model of a House

The Thermostat subsystem models the operation of a thermostat, determining when the
heating system is turned on and off. It can be seen that it is modeled by a single discon-
tinuous Relay block, and outputs either a numerical one or zero, based on the temperature

58
 9 Simulink

error. It is configured to keep the temperature bounded within a specified range from the
set point.
Next, let’s turn our attention to the Heater subsystem which contains some of the most
commonly used blocks. The circular block is a Sum block, which outputs the signed sum
of the input signals. The triangular block is a Gain block, which outputs the product of the
input and the constant defined in the block. The square block with a cross (×) outputs the
product of the two input signals at each instance of time. In the case where one signal is a
constant, it is equivalent to the Gain block.
Interpreting the blocks as an equation, we see that the Heater is modeled as a discontinuous
but static (memoryless) system. The instantaneous heat transfer rate is given by the product
¤ heat capacity 𝑐, and the on/off signal from
of the temperature difference, air flow rate 𝑀,
the Thermostat;

(
𝑀¤ 𝑐 (𝑇Heater − 𝑇room ) , when on
𝑄¤ Heater =
0, when off.

A small comment on the static property of the subsystem — it is assumed that the heating
element reaches the temperature 𝑇Heater immediately when turned on (regardless of the
ambient temperature) and that all heat transfer stops immediately when turned off. This
approximation makes sense when considering that the time unit used in the model is hours
and that the time constant of a conventional heater is likely to be orders of magnitude
smaller — minutes or even seconds.
Finally, the House is modeled as a dynamical system, meaning that the output is determined
not only by the instantaneous input but the entire history of inputs, as well as the initial
conditions. This is achieved using the square Integrator ( 1𝑠 ) block, where 𝑠 is the Laplace
variable. The Integrator is the backbone for building differential equations using block
diagrams, which is how all dynamic systems are modeled in Simulink.
Reading the block diagram as a differential equation, we see that the rate of increase in
temperature of a mass 𝑀 with heat capacity 𝑐 is determined by the heat from the Heater,
and losses – defined by the quotient of the temperature difference inside and outside, and
the equivalent thermal resistance for the whole house 𝑅 eq ;

𝑇room − 𝑇out
𝑀 𝑐 𝑇¤room = 𝑄¤ Heater − .
𝑅 eq

As an optional exercise, you are encouraged to derive this equation, on your own, from the
block diagram representation in Figure 9.6c.

59
 9 Simulink

9.4.3 Simulation resuls

To start the simulation, click the Run button in the Simulation toolstrip, or press F5 . You
can view the results in the Scope block labeled PlotResults (see Figure 9.7). You will see
that the Scope has been partitioned to display the two inputs in separate subplots. The
cumulative heating cost (calculated as the integral of the Heater’s output, multiplied by the
price) should appear in the above plot, and the indoor and outdoor temperatures are shown
in the bottom plot as yellow and blue signals, respectively. By enabling View Legend , you
should see the signals’ labels as well.
You can see that the room temperature oscillates around the set point, determined by how
the thermostat is configured. Note that the cumulative cost remains constant when the
heater is off. Additionally, you can see that the room cools off much faster at night, when
the outside temperature is lower.

Figure 9.7: Simulation results from the Thermal Model of a House

Exercise
Here are a few things to try to see how the model responds to different parameters.
• Open the Constant block labeled Set Point and set the value to some other value.
Rerun the simulation to see how the indoor temperature and heating costs change.
• Adjust the outside temperature (the Avg Outdoor Temp block) and rerun the simula-
tion to see how it affects the indoor temperature.

60
 9 Simulink

• Adjust the daily temperature variation by opening the Sine Wave block labeled Daily
Temp Variation and changing the Amplitude parameter, then re-run the simulation.
• Restore the model to its original state. Then replace the Relay block in the Thermostat
subsystem with a Saturation block, with limits (1,0). The thermostat control signal is
now proportional to the error, and can now assume values between zero (off) and one
(full power). How does this affect the temperature of the house?
• If you zoom in on the curve of the indoor temperature, you may see some small
high-frequency oscillations superimposed on the output. What phenomenon do you
think this is? Try to change the model configuration to remove the oscillations.
Hint: There are nonlinear dynamics! Try a solver which is better at DAEs, see Table 9.1.

9.5 Creating your first dynamic model

In this section, you will learn how to create your own (linear) dynamic system model using
the basic building blocks in Simulink. You should start by closing any models still open
from the previous sections. Then open Simulink, and on the Simulink start page, choose a
new “Blank Model”.
When drawing a block diagram from a differential equation, there are a few rules of thumb
that will help you create a model in a systematic manner.
1. The number of equations is equal to the number of summation blocks you will need.
2. The relative degree (difference between the highest and lowest order of differential
terms) is equal to the number of integrators you will need.
3. Factor out the highest-order term, then start building your model at the Sum block.
For this tutorial, we will consider a linear damped oscillator written as

𝐹 = 𝑚 𝑥¥ + 𝑑 𝑥¤ + 𝑘 𝑥,

where the force 𝐹 is the input to the system, and the deflection 𝑥 is the output.

9.5.1 Walkthrough
Using the rules of thumb listed above, we can identify that we will need one sum block and
two integrators. The next step is to factor the highest-order term. We get

𝑚 𝑥¥ = 𝐹 − 𝑑 𝑥¤ − 𝑘 𝑥,

which tells us that three (not yet present) forces are acting on the oscillator. To get the
acceleration, we need to divide away the mass.
Let’s begin by expressing this using a sum block and a gain block, respectively. Open the
Library Browser, navigate to the Commonly Used Blocks, then drag and drop a Sum and a
Gain block into your blank model.
We can ignore the settings of the Sum block for now, but you should change the Gain value
to 1/m. If the block lights up red, don’t worry about that yet – it is simply a warning that the

61
 9 Simulink

variable m is undefined. Finally, connect the output of the Sum with the input of the Gain
and double-click on the newly created connection to label the signal something appropriate
like mx”. You should have something that looks like Figure 9.8.

mx''

Figure 9.8: Starting at the equals sign

The output signal from the Gain is now the acceleration of the oscillator. To get the
velocity, we should use an Integrator. Instead of using the Library Browser, we will use a
different method. Single-click anywhere in the empty space of the model and begin typing
“integrator”. You should see an auto-completion list appear where you clicked. You can
click the name of the block you want (or hit if it appears at the top of the list) to insert
it. This is the fastest way to insert blocks that you know the name of. Once you have added
the integrator, you can copy ( Ctrl + C ) and paste ( Ctrl + V ) it to get another one. Connect
them in series with the Gain and label the intermediate signals x” and x’, like in Figure 9.9.

mx'' x'' 1 x' 1

s s

Figure 9.9: Adding integrators

Now that we have defined the output signal, we can connect it to a sink. For this example,
visualizing the signal with a Scope will be enough. A Scope can be found under Sinks in
the Library Browser; connect it to the output of the second Integrator, and label the signal
x, see Figure 9.10.

mx'' x'' 1 x' 1 x

s s

Figure 9.10: Adding a Scope

The block diagram is not yet finished. Next, we need to feed back the spring and damper
forces, and then add the input force 𝐹. Let’s turn our attention to the Sum block. By
double-clicking it you can edit the List of signs. The signs are defined counter-clockwise
from the top and can be either plus (+), minus (-), or blank (|). Change it to “-+-”.
Now let’s feed the spring force back onto the mass. We know that the spring force is the
product of the stiffness 𝑘 and the deflection 𝑥, so we need another Gain block. To make
it easy, right-click and hold the 1/m Gain, then drag it somewhere to create a new copy,
Gain1. Change the gain value to k. Then right-click on it and choose Rotate & Flip Flip Block
(or use Ctrl + I ). Place it above the old Gain block and connect the input of the Sum with
the output of Gain1, then connect the input of Gain1 with the signal labeled x. If done
correctly, a branch should be created for you. Rename the signal to kx. You should have
something like in Figure 9.11.

62
 9 Simulink

kx

mx'' x'' 1 x' 1 x

s s

Figure 9.11: Spring force feedback

For the damping force, we will do something much faster. Right-click and hold the signal
labeled x’ to branch a new arrow, then drag it directly to the lower input of the Sum block.
Now right-click and drag one of the old Gain blocks so that the new copy ends up on top
of the newly created line. If done correctly, it should connect itself so that the signal passes
through it. Change the gain value to d and rename the signal, see Figure 9.12.

kx

mx'' x'' 1 x' 1 x

s s

dx'

Figure 9.12: Damper force feedback

Finally, let’s add an input signal. You can find a Step block under the Sources tab in the
Library Browser. Set the Step time to 0, and the Final value to F. The final block diagram
can be seen in Figure 9.13.

kx

F mx'' x'' 1 x' 1 x

s s

dx'

Figure 9.13: Adding an input signal

The last thing to do before we can run the model is to define the parameters. There are a
few ways to do this, the simplest way is to define them through the Command Window so
that they are temporarily saved in the Workspace.

63
 9 Simulink

>> F = 1; % Force [N]

>> m = 1; % Mass [kg]
>> k = 1; % Stiffness [N/m]
>> d = 1; % Damping [Ns/m]

Finally, you can test the simulation by running it. The output should look like Figure 9.14.

Figure 9.14: Step response of the damped oscillator model

If the model appears correct, you can save it as [Link] – it will be used for
the following sections.

9.5.2 Partitioning the Scope

Oftentimes, it is convenient to visualize multiple signals inside the same scope, to see how
they interact during the evolution of the system. A fundamental rule that should not be
broken is: do not mix signals of different units in the same axes! Instead, you can partition
the Scope window, similar to the subplot command. A good idea is to partition it along
the vertical so that the signals in the subplots share the same time axis. Have a look at
Figure 9.7, which does exactly this in order to separate cost from temperature.
Now turn your attention to the damped oscillator model. Click and hold the right mouse
button at the corner of the signal labeled x’ to create a new branch, and connect it to
the scope. If done correctly, Simulink will automatically increase the number of ports to
facilitate the connection; see Figure 9.15.

64
 9 Simulink

kx

x
F mx'' x'' 1 x' 1
s s

dx'

Figure 9.15: Adding velocity to the Scope

Finally, open the scope, click the cogwheel in the toolbar to bring up the Configuration
Properties, and click Layout. Choose the 2 × 1 grid configuration (Figure 9.16), then accept
the changes and close the configuration window.

Figure 9.16: Configuration to partition the Scope

After running the simulation again, the Scope should look like Figure 9.17.
For more practice on building a block diagram based on a system of differential equations,
see this tutorial from the Control Tutorials for MATLAB and Simulink.

65
 9 Simulink

Figure 9.17: Step response of the damped oscillator, showing both deflection and velocity

66
 9 Simulink

9.6 Matlab integration

As shown in Section 9.4.2, Simulink models are hierarchical – where each subsystem’s
inputs and outputs inhabit the level above it. When working in the first layer of a model,
the layer “above” is the Matlab development environment and its Workspace. This section
will show how to leverage this connection to treat Simulink models more like Matlab
functions, with a high level of control over the model parameters and output signals. The
first thing that you will learn to control is how run-time data can be exported from Simulink
to the Workspace.

9.6.1 Exporting outputs

By default, Simulink will export the simulation time vector, as well as any outputs defined
in the topmost level of the model hierarchy. The data can be exported in various formats,
we will use the simplest, which is just a numerical array. Since Matlab R2019a, the default
behavior is also to group the exported values into a single Simulation Output data structure
typically named out, we will disable this feature for now.
You can configure the data to export through the Model Settings (you can get there from
the Modeling Toolstrip, under the Setup tab) then navigate to the Data Import/Export tab.
For the sake of this tutorial, make sure that [Link] has the same settings
as shown in Figure 9.18.

Figure 9.18: Data export configuration for the damped oscillator model

67
 9 Simulink

Subsequently, we need to add some outputs. Add two Out1 (sink) blocks, using any method
you prefer, and connect them to signals x and x’, so that the deflection has number 1 and
the velocity has number 2, see Figure 9.19.

kx
1

x
F mx'' x'' 1 x' 1
s s

dx'
2

Figure 9.19: Adding outputs

Finally, run the simulation again and wait until it finishes. Check the Matlab Workspace.
There should be two vectors named tout and yout, which are the time and output vectors,
respectively. Each column in yout corresponds to an Out1 block, with the ordering given
by the respective port numbers.
You can try the following commands to plot the signals using Matlab.
>> subplot (2 ,1 ,1) , plot(tout , yout (: ,1)) % Plot deflection
>> subplot (2 ,1 ,2) , plot(tout , yout (: ,2)) % Plot velocity

The resulting figure should be equivalent to the output in the Scope.

9.6.2 To Workspace
Furthermore, there is another sink called To Workspace, which allows you to export in-
dividual signals as various data types. This allows greater flexibility when you have a
large model and may want to export internal signals, such as controller logic, which aren’t
considered model outputs. In the case where the Single simulation output option is on, the
saved variables will appear inside the output structure.

9.6.3 Running simulations

Now that we have seen how to configure which model outputs, we can look at how to call
model files from the Command Window, similar to functions. The command to use is sim
and is called with the path to the model as a string. If you do not provide the full path, you
must be working in the correct directory
>> sim('DampedOscillator ')

This command is equivalent to clicking the Run button of the model, and will save the out-
puts in the Workspace as specified in the Configuration Parameters previously. However,
it is also possible to specify the variables in which to save the outputs.

68
 9 Simulink

Note that when using this syntax, the returned data is not determined by the data export
configuration, and you have to make sure that you are saving the correct signals!
For example, you can save the time vector, internal states, and output matrix as t, x and y
respectively, using
>> [t,x,y] = sim('DampedOscillator ') % Run Simulink model

9.6.4 Single simulation output

When working with larger models, it is recommended to use the single output option, to
prevent erroneous variable labeling. Open the data export configuration of the model, and
tick the “Single simulation output” option. Then try the following command
>> simout = sim('DampedOscillator ') % Run with single output

You will need to use dot-notation to access the exported data, and keep in mind that like
most things in Matlab, the output structure can be indexed. For a review of structures,
see Section 3.6.

Exercise
In this exercise, you will try your hand at running the damped oscillator model with varying
model parameters, and plotting the results within the Matlab IDE.
Create a new Matlab script, and add the following model parameter initialization
clear all , close all , clc

% Constant parameters
F = 1; % Force
m = 1; % Mass
k = 1; % Stiffness
d = 1; % Damping

Then add the following functionality:

1. Add the command to run the damped oscillator model, saving the output structure.
2. Plot the deflection in a Figure in Matlab, and make sure it is correct.
3. Now, change the code so that the model is simulated four times, each time with a
different damping value. You can try 𝑑vec = [2.0, 1.0, 0.5, 0.2] for instance.
Hint: You may want to consider using iterations, see Section 4.3.
4. Plot the results in four separate Figures.
5. Label the plots accordingly—-so that it becomes clear which damping values were
used.
Hint: You can read more about formatting text on the MathWorks homepage.
6. Optional: Export the Figures as separate vector graphics files, see Section Section 6.4.

69
 9 Simulink

9.7 Dynamic system models

At the beginning of this Chapter, we reviewed some of the different sub-headings of the
Simulink Library Browser. In the Continuous and Discrete libraries, you will find blocks
for specifying dynamic system model objects (State-Space, Transfer Fcn, PID-Controller,
etc.) directly in Simulink, in continuous and discrete time, respectively. The parameters for
these blocks take on a form very similar to that of their corresponding Matlab commands
in the Control System Toolbox (presented in Chapter 7).

9.7.1 The Transfer Fcn block

For the sake of comparison, let’s examine how to simulate the damped oscillator model as
a transfer function instead. If you were instantiating the model in Matlab, you would have
used the tf-command – the procedure in Simulink is similar.
It won’t be shown here, but the transfer function corresponding to the differential equation
from Section 9.5 is given by

𝑋(𝑠) 1
𝐺(𝑠) = = .
𝐹(𝑠) 𝑚𝑠 +𝑑𝑠+𝑘
2

Open the damped oscillator model, and add a Transfer Fcn block to the model, using some
method you prefer. Upon adding the block, open the parameters and specify the following
• Numerator coefficients: [1]
• Denominator coefficients: [m, d, k]
Next, add a new Scope, and make sure it has 2 inputs and is partitioned like in Section 9.5.2
– if you create a copy of the previous Scope, the settings will carry over. Finally, connect
it to both the output of the Transfer Fcn and the deflection signal x, and connect the input
of the Transfer Fcn to the Step output by creating a new branch. Your model should look
something like Figure 9.20.

kx
1

x
F mx'' x'' 1 x' 1
s s x

dx'
2

Figure 9.20: Comparing the output of a Transfer Fcn to the block diagram

70
 9 Simulink

Upon running the simulation, you should be able to convince yourself that the two models
produce “the same” (sufficiently similar) outputs. As an optional exercise, modify the
Simulink model so that the new Scope displays the difference between the two outputs in
a third subplot. What is the order of magnitude of the difference?

9.8 Closing words

Hopefully, by now, you grasp enough of the fundamentals that using new commands,
blocks and Toolboxes correctly now becomes a matter of underlying theoretical knowledge
and not a question of understanding the basic framework that is Matlab/Simulink.
The author of this text hopes you enjoyed working through this tutorial (or that you at
least found it useful) and wishes you the best of luck with your future endeavors with the
modeling and simulation of dynamic systems using Matlab/Simulink.

71
Solutions to exercises

Elementary mathematics
clear all , close all , clc

% Variable x
x = pi;

% First set of functions

f1 = cos(x) + i*sin(x);
f2 = exp(i*x);

% Second set of functions

g1 = 1/(1+ x);
g2 = 0.5 / (1+i*sqrt(x)) + 0.5 / (1-i*sqrt(x));

% Test equality
diff_f = f1 - f2
diff_g = g1 - g2

Comments
Even though the functions ( 𝑓1 (𝑥), 𝑓2 (𝑥)) and (𝑔1 (𝑥), 𝑔2 (𝑥)) are in fact equivalent, and thus
should be equal for all values 𝑥 ∈ ℂ, there are some values for which the difference is
non-zero. One such case is for 𝑥 → −1, for which 𝑔1,2 (𝑥) → ∞. Matlab solves the equation
as
Inf − Inf = NaN (not a number).

Additionally, you will encounter cases where the answer is in the order of magnitude
1×10=16 . This is a numerical residual, and in most cases can be considered zero.

72
 Solutions to exercises

Vectors and matrices

clear all , close all , clc

% Declare vector and number of rows

a = [1, 2, 3];
m = 3;

% Number of columns
r = size(a ,2) - 1;
n = m + r;

% Diagonal columns
D = repmat (a, m, 1);

% Index vector
R = 0:r;

% Sparse diagonal
A = spdiags (D, R, m, n);
full(A)

73
 Solutions to exercises

The Control System Toolbox

clear 'all ', close 'all ', clc

% Two ways to declare the transfer functions

s = tf('s');
G = 100/( s^2 + 10*s + 100)
C = tf ([4 ,40 ,400] ,[1 ,40 ,0])

% Closed -loop system

G_cl = feedback (G*C ,1 , -1)

% Poles and zeros

zpk(G_cl) % The term (s^2 + 10s + 100) can be cancelled

% Minimal realisation
G_cl = minreal (G_cl);

% PID controller
C = pid(C)

% Linear analysis
damp(G), stepinfo (G)
damp(G_cl), stepinfo (G_cl)

% Figures
figure ; step(G,G_cl);
legend ('Uncontrolled ','Closed -loop ')
print('ExerciseDynamicSystemStep .eps ','-depsc ','-painters ','-loose
')

figure ; pzmap (G,G_cl);

print('ExerciseDynamicSystemPoleZero .eps ','-depsc ','-painters ','-
loose ')

Comments
In the uncontrolled system, the damping ratio of 0.5 gives rise to oscillations, with a
maximum overshoot of 16.29 %. In the controlled system, the two poles are moved onto
the real axis, which corresponds to a damping of 1, and so no oscillations or overshoot is
present. This, paired with the fact that the poles lie further to the left, gives a faster settling
time. The resulting plots can be seen in Figure 1.

74
 Solutions to exercises

Step Response
Pole-Zero Map
1.2
10
Uncontrolled
Closed-loop 8
1
6
Imaginary Axis (seconds-1 )

0.8 4
Amplitude

2
0.6
0

-2
0.4
-4

0.2 -6

-8
0
0 0.2 0.4 0.6 0.8 1 1.2 -10
-25 -20 -15 -10 -5 0
Time (seconds)
Real Axis (seconds -1 )

(a) Step response (b) Pole-zero map

Figure 1: Linear analysis plots of the controlled and uncontrolled system

75
 Solutions to exercises

Discrete-time systems
clear 'all ', close 'all ', clc

% State space matrices

A = [0, 1; -1, -1];
B = [0; 1];
C = [1, 0];

% Continuous -time system model , for comparison

sysC = ss(A,B,C ,0);

% Appropriately sized identity matrix

I = eye(size(A));

% Repeat the simulation for different orders of approximation ...

for N = [1, 2, Inf]

% ... and different sample times

for Ts = [0.5 , 0.1 , 1.0]

% Calculate the discrete matrices for the different

approximations
switch N
case 1, Phi = I + A*Ts;
case 2, Phi = I + A*Ts + 1/2*( A*Ts)^2;
otherwise , Phi = expm(A*Ts);
end
Gamma = A\(Phi -I)*B;

% Approximated system model

sysA = ss(Phi ,Gamma ,C,0,Ts);

% Discretised system model , using ZOH

sysD = c2d(sysC ,Ts ,'zoh ');

% Comparison
figure ('Name ',sprintf ('N = %d, Ts = %1.1f', N, Ts));
step(sysC ,sysD ,sysA);
legend ('Continuous ','Discretisation ','Approximation ')
print ( sprintf ('ExerciseDiscretisation -N%d-T%[Link] ',N,Ts),
'-depsc ','-painters ','-loose ');
end
end

76
 Solutions to exercises

Comments
As can be seen in Figure 2, the “accuracy” of the approximation is generally more important
than the step size. For the first-order approximation, choosing a sample-time of 1 sec
will make the system unstable, whereas using the same sample time in the second-order
approximation captures the dynamics quite well.
We can identify the approximation as the truncated Taylor expansion for the exponential
function
1 1
exp(𝑥) = 1 + 𝑥 + 𝑥 2 + 𝑥 3 + . . .
2! 3!
Using the matrix exponential command expm, we see that ZOH is equivalent to taking the
matrix exponential of the system matrix multiplied with the sample time.

Step Response
Step Response
1.5
1.2
Continuous
Continuous
Discretisation
Discretisation
Approximation
1 Approximation

1
0.8
Amplitude

Amplitude

0.6

0.5
0.4

0.2

0
0 5 10 15 20 25 0
0 2 4 6 8 10 12 14 16 18 20
Time (seconds)
Time (seconds)

(a) 𝑁 = 1, 𝑇s = 0.5 sec (b) 𝑁 = 2, 𝑇s = 1.0 sec

Figure 2: Comparison between order of approximation and step size

77
 Solutions to exercises

Simulink-Matlab integration
clear 'all ', close 'all ', clc

% Constant parameters
F = 1; % Force
m = 1; % Mass
k = 1; % Stiffness

% Define damping vector to iterate through

d_vec = [2.0 ,1.0 ,0.5 ,0.2];

for n = 1: length ( d_vec )

% Choose damping value

d = d_vec (n);

% Run Simulink model , enforcing the single - output struct

out(n) = sim('DampedOscillator ','ReturnWorkspaceOutputs ','on ');

% Plot results
figure (n);
plot(out(n).tout ,out(n).yout (: ,1))
grid on; set(gca ,'YLim ' ,[0 ,2]);

% Add text , print Figure

title( sprintf ('Step response with d = %1.1f',d));
xlabel ('Time ( seconds )'); ylabel ('Amplitude ');
print( sprintf ('ExerciseSimulink -%1.1 [Link] ',d), ...
'-depsc ','-painters ','-loose ');
end

Comments
We can see that the system becomes critically damped at 𝑑 = 2, and decreasing the damping
further introduces oscillations in the system, see Figure 3. It follows that a system with no
damping (𝑑 = 0) will not display any attenuation in the amplitude of the oscillations.

78
 Solutions to exercises

Step response with d = 2.0

2 Step response with d = 1.0
2
1.8
1.8
1.6
1.6
1.4
1.4
1.2
Amplitude

1.2

Amplitude
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0 1 2 3 4 5 6 7 8 9 10 0
0 1 2 3 4 5 6 7 8 9 10
Time (seconds)
Time (seconds)

(a) 𝑑 = 2.0 (b) 𝑑 = 1.0

Step response with d = 0.5 Step response with d = 0.2
2 2

1.8 1.8

1.6 1.6

1.4 1.4

1.2 1.2
Amplitude

Amplitude

1 1

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

0 0
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
Time (seconds) Time (seconds)

(c) 𝑑 = 0.5 (d) 𝑑 = 0.2

Figure 3: Running the damped oscillator model with different damping values

79
Appendix

Creating a publication-quality Figure

%% Preamble
clear 'all ', close 'all ', clc % Clean slate

% Custom colors
c1 = [31 , 119 , 180]/255; % #1 f77b4
c2 = [255 , 127 , 14]/255; % # ff7f0e

% Data to plot
T = 1; % Time constant (sec)
t0 = 0; tf = 6*T;
t = linspace (t0 ,tf ,50); % Time vector for simulation
f1 = 1 - exp(-t./T); % Exponential smoothing
f2 = exp(-t./T)/T; % Exponential decay
tT = [t0 ,t0+T]; yT = [0 ,1]; % Help - lines for gradient at t=0

%% Figure dimensions and placement

% Figure size in cm
width = 12;
height = width * 2/(1+ sqrt (5));
pos = [0, 0, width , height ];

% Screen placement
set(groot ,'Units ','centimeters ')
scrn = get(groot ,'ScreenSize ') - pos;
pos (1:2) = scrn (3:4) /2;

%% Create Figure programmatically

% Add Figure window
fig = figure ('Units ','centimeters ','Position ',pos);

% Add Axes
ax = axes('Parent ',fig ,...
'XGrid ','on ','YGrid ','on ','XColor ','k','YColor ','k',...

80
 Appendix

'XMinorTick ','on ','YMinorTick ','on ','YLim ' ,[0 ,1.2] ,...

'TickLabelInterpreter ','latex ','FontSize ' ,10);

% Set axis labels

set([Link] ,'String ','Time (sec)','Interpreter ','latex ');
set([Link] ,'String ','Amplitude ','Interpreter ','latex ');

% Add Lines
ln (1) = line(ax ,t,f1 ,'Color ',c1 ,'LineWidth ' ,1.0,'LineStyle ','-');
ln (2) = line(ax ,t,f2 ,'Color ',c2 ,'LineWidth ' ,1.0,'LineStyle ','-');
ln (3) = line(ax ,tT ,yT ,'Color ','k','LineStyle ','--');

% Add Legend
leg = legend ('Parent ',fig ,...
'String ',{'$f(x) = 1-e^{-t/T}$','$f''(x) = 1/T\cdot e^{-t/T}$'},
...
'Interpreter ','latex ','FontSize ' ,10,'Location ','southeast ');

%% Print Figure to file

% Resize Axes to fill Figure
p = get(ax ,'TightInset '); % Read only!
set(ax ,'Position ',[p(1) ,p(2) ,1-p(1) -p(3) ,1-p(2) -p(4) ]);

% Print Figure to file , using the vector format

print(fig ,'-depsc ','-painters ','PublicationQualityFigure .eps ');

% Optional : Convert to PDF ( tested with MikTeX )

! epstopdf PublicationQualityFigure .eps --outfile ="
PublicationQualityFigure -eps -converted -[Link]"

81
 Appendix

Custom linear analysis plots

clear 'all ', close 'all ', clc

% Second - order system model

o = 1; d = 0.2;
sys = tf ([o^2] ,[1 ,2*d*o,o^2]);

for i=1:4

% Create a new Figure for each plot

fig(i) = figure (i);

switch i
case 1
% Step response
obj(i) = stepplot (sys);
case 2
% Pole -zero map
obj(i) = pzplot (sys);
setoptions (obj(i), 'Grid ','on ',...
'XLim ' ,[ -1.8 ,0] , 'YLim ' ,1.5*[ -1 ,1]);
case 3
% Bode diagram ( magnitude only)
obj(i) = bodeplot (sys);
setoptions (obj(i), 'Grid ','on ', 'PhaseVisible ','off ');
case 4
% Nyquist diagram
obj(i) = nyquistplot (sys);
setoptions (obj(i), 'Grid ','on ');
end

% Extract the Title from each Plot

opt = obj(i). getoptions ;
s = opt. Title . String ;

% Remove all forbidden characters

s( regexp (s,'[-, ]'))=[];

% Use the plot title as part of the filename

filename = ['LinearAnalysis ',s,'.eps '];

% Print Figure to file , using the vector format

print(fig(i),'-depsc ','-painters ','-loose ',filename );

end

82