AIAA-2000-4395
VECTORIZATION BLOCKS FOR SIMULINK
Giampiero Campa‡, Mario Innocenti†,
Department of Electrical Systems and Automation, University of Pisa, 56126 Pisa, ITALY
Abstract
The paper presents two blocks that extend
Simulink range of applicability to varying matrices,
allowing the user to build (without writing a single
line of code) dynamical systems that can be both
complex and suitable for real time simulations. First
of all, the three common techniques to build
dynamical system with Simulink will be presented,
showing the advantages and disadvantages of each of
them. Then we will explain why (when dealing with
complex multivariable systems) neither of the three
can be at the same time easy to use and suitable for
real time simulations. Two Simulink blocks will then
be proposed that solve this problem by allowing
varying matrices to be handled directly in the
Simulink environment.
The aspect, the user
interface, and the most important part of the code of
these blocks will be discussed, and their use will
finally be explained by means of some examples. In
particular, we will present an example on how to
implement a complex multilayer adaptive neural
network using these two blocks.
Introduction: S-functions
S-functions are the way in which a dynamical
system is represented inside Simulink environment.
Every Simulink object (block, composition of
blocks, schemes) is a dynamical system; hence,
every object (included the entire simulation scheme)
is represented by an S-function.
Typically, the built in Simulink blocks are associated
with precompiled S-functions, while user added Sfunctions are obtained joining blocks or directly
written in a language such as Matlab, C, C++,
Fortran.
______________________________
†Professor, Associate Fellow AIAA
‡Research Assistant Professor, currently with Department of
Mechanical and Aerospace Engineering, West Virginia
University, Morgantown, WV, USA
Copyright © 2000 by Mario Innocenti and Giampiero Campa,
Published by the American Institute of Aeronautics and
Astronautics, Inc., with permission
Matlab and Simulink are trademarks of The Mathworks Inc.
Whatever the language in which the S-function is
written in, it must interact with Simulink in the same
way, that is, it must have the structure shown in [1].
Different ways to build S-functions
Joining together Simulink blocks
Perhaps the most common way to build
dynamical systems (hence S-functions) in Simulink,
is just to join together Simulink blocks. For example,
if we have two blocks that represent two systems and
we would like to have the system consisting of the
series of the two, the easiest way to obtain the series
system is just drawing a wire that connects the
output of the first system to the input of the second
one.
We can then group the two systems together using
the Simulink “Group” function in order to obtain a
block that appears as any Simulink block and
represents our “new” dynamical system.
The ease of operations like the one depicted
above has been maybe the principal cause of
Simulink success as an interface to deal with
dynamical systems, indeed in this way we can build
a dynamical system straightforwardly without the
need to write a single line of code.
In some cases however assembling subsystems
may not be the best choice, this is true especially if
we already have the input-state-output equations of
the system that we are going to build, and especially
for nonlinear multivariable systems.
Let’s consider for example a system in which the
output depends from the state vector in this way:
y = C(x )D(x )x
where each entry of the matrices C and D is a given
function of x. If we want to obtain y by connecting
some blocks that transform the signals x into y, then
the only solution is to do the multiplication in a wire
by wire fashion, if both C and D are 2 by 2 matrices
we have:
y1 = (C11 ( x ) D11 ( x ) + C12 ( x ) D21 ( x )) x1
+ (C11 ( x ) D12 ( x ) + C12 ( x ) D22 ( x )) x2
y2 = (C 21 ( x ) D11 ( x ) + C 22 ( x ) D21 ( x )) x1
+ (C 21 ( x ) D12 ( x ) + C 22 ( x ) D22 ( x )) x2
1
American Institute of Aeronautics and Astronautics
(1)
For instance, if we use the two blocks C(x) and D(x)
to construct the respective matrices, the resulting
Simulink diagram is:
stored in the variables C and D, it takes just the
simple expression
y=C*D*x;
1
x
Selector
[C11 C12]
Dot Product 1
C
Selector
Dot Product
Dot Product 2
[C21 C22]
Dot Product 4
[D11 D21]
Dot Product 5
[D12 D22]
C(x)
1
Selector
y
to compute the output. If we save this file as
“mysys.m” then we can use it just by specifying the
name “mysys” in the mask of the predefined
Simulink S-function block:
D
Selector
Dot Product 3
D(x)
1
y
oursystem
S−Function
1
x
Figure 1. Matrix Multiplication Example
Figure 2. Simulink S-function block
Needless to say, this scheme is not general nor it
is easily adaptable if the dimensions involved vary,
and especially if they rise above 3 or 4. Where does
all this complexity come from?
The main problem in this case is that, although
the wires at the output of the blocks C(x) and D(x)
carry all the four elements of the two matrices, when
we have to multiply the two matrices we have to
decompose them into rows or columns. In
conclusion, even though constructing systems by
connecting blocks has in general several advantages,
in some cases (especially when we have to handle
nonlinear multivariable systems) it gets very
difficult.
Writing S-functions as M-files
Another option is to write directly the S-function
that we need using Matlab language. Since Matlab
notation is very similar to the standard mathematical
notation, this method, when properly used, naturally
leads to the most elegant and compact way to
describe dynamical systems. For example, in the
case of the system considered above, we could write
the following M-file:
function
[y,x0,st,ts]=mysys(t,x,u,flag)
if flag==0
% initial info
% 0 continue and discrete states
% 2 inputs, 2 outputs
y=[0,0,2,2,0,0,1];
x0=[];
st=[];
ts =[0 0];
% sample time
elseif flag == 3 % outputs
C=… {compute C(x)}
D=… {compute C(x)}
y=C*D*x;
else
y=[];
end
Without going into the details we will just point out
that once we have C(x) and D(x) computed and
This block will then be equivalent to the scheme in
Figure 1. Although we no longer “see” how the
system is built directly from the Simulink scheme,
when dealing with complex systems this method
seems to be by far the simplest one. Writing Sfunctions as M-files has only one serious drawback:
since we use Matlab language the Matlab interpreter
has to be called to evaluate each row of the Sfunction at every simulation step, this fact prevents
the scheme from being used with Real Time
Workshop, (which is the standard real time code
generator for Matlab) and, last but not least, may
slow down the simulation from 10 to 50 times.
There exist compilers that translate M-files to C-files
and then to mex (Matlab executable) files, but
unfortunately none of them (yet) produce a code that
can be used with Simulink, so the only way in which
we can use such compilers is the following:
1. write the main S-function as a frame that in
turn calls other subfunctions written as Mfiles themselves.
2. use the compiler on these subfunctions.
This method can help to speed up the simulation
when the subfunctions carry out a significant amount
of computation, but since we still need to call the
interpreter, we can’t use Real Time Workshop yet.
Writing S-functions as C-files
Instead of using Matlab language, we could write
the S-function in C language directly, following the
“Simulink level 2” conventions [1]. In this way, the
S-function could be compiled so that it behaves as a
built-in Simulink block, so the interpreter no longer
needs to be called, Real Time Workshop can be
used, and the simulation is as fast as possible.
Unfortunately, writing an S-function as a C file is
in no way as simple as writing it as an M-file, and
requires, other that a good C programming
knowledge, a careful study on how Simulink works,
on how to pass variables and parameters from the
Simulink block to the inner s-function, and so on.
But the biggest problem is that the Matlab language,
2
American Institute of Aeronautics and Astronautics
and its types are no longer available so we cannot
rely on indexing, built in functions, and above all we
cannot rely on matrix operations.
In the above example, the only way for us to
perform the output computation with a line of code
as simple as it was in the Matlab case, i.e.:
y = C*D*x;
would be to define a class matrix, a subclass vector,
and overload the operator *, otherwise we have to
write a function mat_product that performs the
matrices multiplication, and write something like:
mat_product(tmp1,D,x,2,1);
mat_product(y,C,tmp1,2,1);
the same is true if we want to perform in C any other
Matlab built in operation.
In conclusion, writing the S-function in C it is often
the best thing from the performance point of view
but unfortunately it’s not very practical and requires,
some programming knowledge and a considerable
amount of work.
Vectorization Blocks
Main Point
The signals that flow through Simulink wires,
which in general are the only quantities allowed to
vary with time, can be scalars or vectors.
It is not permitted to carry a matrix through a
wire, namely we could certainly carry a p by n
matrix as a vector of p*n entries (as we do in Figure
1 with the matrices C(x) and D(x)) but all the built in
blocks that will process this vector will still treat it
just like a vector, so every time we have to process
the p by n matrix we have to decompose it into its
elements. These decompositions, and the successive
compositions make it very uncomfortable to deal
with varying matrices in the Simulink environment.
All this could be avoided if we had some blocks that
would be capable to deal with matrices without
decomposing them into elements.
This is exactly the main idea behind the
“vectorization blocks”. Each input of them is a
vector that represents a matrix ordered columnwise
according to dimensions provided through its mask,
(so row vectors, column vectors and even scalars are
just a particular case), and the output is another
vector that represents the “output matrix”, which is
the result of some operations on the input matrices,
(product, transposition, norm, pseudo-inversion and
so on).
These blocks are based on S-functions written in C
(following the Simulink 3 optimized “level 2”
conventions) and then compiled so that their
execution is as fast as if they were built in blocks.
We will present the blocks for matrix transposition
and matrix multiplication, since they can cover the
implementation of the majority of the control
systems. The blocks for pseudo inversion, inversion,
norm, determinant, (the inversion block would very
useful for simulating non-rigid bodies), are currently
under construction, that is, these blocks exists but
they are not yet usable within real-time workshop.
Multiplication Block
The first block, which is based on the S-function
“vmult.c”, performs matrix multiplication. This
block accepts as inputs two vectors each one
representing a matrix whose elements are taken
columnwise, and accepts as parameters 3 numbers:
1. n = number of rows of the first matrix.
2. p = number of columns of the first matrix =
number of rows of the second one.
3. m = number of columns of the second
matrix.
The two matrices are then reconstructed from the
two vectors, and the first one is (left) multiplied by
the second one, the elements of the n by m resulting
matrix are then reordered columnwise into the output
vector. The following picture shows the appearance
of the block:
x
1
Z=X*Y
2
y
1
z
Vmult
Figure 3. Multiplication Block
By double clicking on the block the following mask
appears:
Figure 4. Multiplication Block Mask
The vector [n p m] contains the dimensions that the
function needs to build the matrices from the inputs.
Transposition Block
The second block, which is based on the Sfunction “vtrsp.c”, performs matrix transposition.
This block accepts as input one vector representing a
3
American Institute of Aeronautics and Astronautics
matrix whose elements are taken columnwise, and
accepts as parameters 2 numbers:
1. n = number of rows of the matrix.
2. m = number of columns of the matrix.
The matrix is then reconstructed from the two
vectors, and then the elements of its transpose are
reordered columnwise into the output vector. The
following picture shows the appearance of the block:
1
Y=X’
1
x
vtrsp
y
real_T *y
= ssGetOutputPortRealSignal(S,0);
InputRealPtrsType up1
=
ssGetInputPortRealSignalPtrs(S,0);
InputRealPtrsType up2
=
ssGetInputPortRealSignalPtrs(S,1);
for(i = 0; i < n[0]; i++)
for(j = 0; j < n[2]; j++)
for(y[i+j*n[0]]=0,k=0;k<n[1];k++)
Figure 5. Transposition Block
y[i+j*n[0]]+=(*up1[i+k*n[0]])*
By double clicking on the block the following mask
appears:
Figure 6. Transposition Block Mask
The vector [n m] contains the dimensions that the
function needs to build the matrix from the input.
Redundancy
Knowing the number of inputs, the number of
dimensions to be typed in the masks of both blocks is
redundant by 1. We choose to retain the present level
of redundancy because it allows double-checking
and therefore helps preventing errors in connecting
the blocks to one another.
Implementation Code
Without going into the details on how to write an
S-function in C, we think it could be useful to have a
quick look at the heart of the two S-functions: the
sub-function “mdlOutputs”, which computes the
output from the values of inputs.
Vmult.c
Here is the code of mdlOutput for the S-function
vmult.c:
static void
mdlOutputs(SimStruct
tid)
{
int_T
i, j, k, *n
=ssGetIWork(S);
*S,
int_T
(*up2[k+j*n[1]]);
}
A quick explanation of the code is necessary.
Integers i, j, k, are indexes, n is a pointer at the three
elements vector of the dimensions, y is a pointer at
the output vector, up1 and up2 are pointers at the
input vectors. The three for loops perform a matrix
multiplication, in particular the first two loops select
the row and column of each element of the output
matrix, and the inner loop performs the scalar
product between the relative row of the left matrix
and the relative column of the right matrix.
It is important to notice that we are deliberately
using a fast indexing that allows us to handle
matrices directly in their vectorized form without
actually losing computation time on array
construction. Thus the code is very optimized versus
speed.
Vtrsp.c
Here is the code of mdlOutput for the S-function
vtrsp.c:
static void
mdlOutputs(SimStruct
*S,
int_T
tid)
{
int_T i, j, k, *n
= ssGetIWork(S);
real_T *y
= ssGetOutputPortRealSignal(S,0);
InputRealPtrsType up
=
ssGetInputPortRealSignalPtrs(S,0);
for(i=0;i<n[0];i++)
for(j=0;j<n[1];j++)
y[j+i*n[1]]=(*up[i+j*n[0]]);
}
Integers i, j, k, are indexes, n is a pointer at the
two elements vector of dimensions; y is a pointer at
4
American Institute of Aeronautics and Astronautics
the output vector, up points to the input vector. The
two for loops perform the matrix transposition just
by exchanging indexes. As in the preceding case,
matrices are handled directly in their vectorized
form.
output layer) and N (connecting input layer to hidden
layer). The output of the hidden layer is the vector
Examples and Considerations
where n1 … nh are the rows of NT and σ is the
activation function, which is chosen to be:
Two multiplications
Let us reconsider the case in which we have the
output depending from the state vector in this way:
σ ( N T x ) = [1 σ 1 ( n1 x1 ) ... σ h ( n h x h )]T
(2)
1
1 + e −αz
σ i ( z) =
(3)
y = C(x)*D(x)*x
we can construct this function using the Varying
Matrices Multiplication Block in this simple way:
C
Z=X*Y
1
Z=X*Y
y
C(x)
vmult 1
vmult 2
D
D(x)
1
x
Figure 7. Two matrix Multiplication with vmult
A quick comparison between Figure 1 and Figure
7 is sufficient to appreciate the simplification;
moreover, this new Simulink scheme is completely
general, that is, if we change the dimensions of the
matrices involved, we only have to write the new
dimensions into the two Varying Matrices
Multiplication Blocks.
The matrix σ’(NTx) is the Jacobian of σ(NTx) with
respect to NTx evaluated at NTx. Z is defined as
diag(M,N), and ZF is its Frobenius norm. All other
values are constants. Specifically F and G are the
learning rates of the first and second layers, and λ is
commonly known as the e-modification gain.
Constants required for the robustifying term are Kv
and Kz. Finally, Z is an upper bound for ZF.
The following picture, Figure 8, shows
realization of the net as a Simulink scheme.
[N]
Y=X’
[x]
N’
[1 0]’
Z=X*Y
[Ntx]
bias
Switch
a
a
[Ntx]
m
[sig]
m
[sigp]
a
1
asig(1−sig)
[sigp]
[Ntx]
sigpNtx
− G[ x ζM T σ ' ( N T x ) + λ ζ N ]
N&
&
T
T
T
−
−
+
F
[(
(
N
x
)
'
(
N
x
)
N
x
)
M
]
σ
σ
ζ
λ
ζ
M
=
T
T
vad
M σ (N x)
K v ( Z F + Z + K Z )ζ
v
m
N’x
Z=X*Y
[sig]
[z]
MIMO Adaptive Sigmoidal Neural Network
The following example is more complex, here we
used the two blocks to implement a continuous time
adaptive multilayer MIMO sigmoidal neural network
[2], [3], [4] which is represented by the following
equation:
the
(sig−sigp)z
−F
[M]
1/s
4
[M]
M
lam n(z) M
Dot Prod
[z]
lam
sqrt
[N]
−G
lam n(z) N
[x]
[M]
Y=X’
[z]
2
[x]
Z=X*Y
zM’
M’
z
1
[N]
1/s
3
N
Z=X*Y
xzM’sigp
[sigp]
zM’sigp
x
[sig]
Z=X*Y
M’sig
M’
1
vad
[z]
Y=X’
[M]
Z=X*Y
M
N’
[N]
M’M
K
trace(M’M)
sqrt
Y=X’
Frob norm
Kv
Z_bar
Z=X*Y
N
N’N
2
K
trace(N’N)
Kz
Product
vbar
Figure 8. Neural Network Scheme
There are two inputs to the net, the first is x,
which should contain information necessary to the
net to approximate a function ∆(x), and the second is
ζ, which is the “error” term that actually drives the
network.
The network output is composed of two vectors:
the adaptive term ν ad and the robustifying term ν ,
which is made proportional to ζ. The network states
are exactly the interconnection weights between
different layers, namely all the entries in the two
synaptic matrices M (connecting hidden layer to
The upper part of the scheme is devoted to the
computation of the hidden layer output (4) and its
jacobian matrix. The central part is the heart of the
function in that it computes the derivatives of the
states N and M, and then integrates them. Finally the
lower part of the scheme contains the blocks for the
output computation, which involves also the
computation of the Frobenius norm of Z.
This scheme has been successfully used to
simulate in real time several neural network based
5
American Institute of Aeronautics and Astronautics
controls [2], [3], although slightly faster, its version
as a C program was 420 lines long and had several
limitations, a full-general and optimized version
would probably be much longer.
Rotation Matrices Library
The multiplication and transposition blocks were
used to build a library of rotation matrices, based
both on Euler angles and quaternions. These matrices
are a classical example of very required varying
matrices blocks, since they are very useful for the
simulation of rigid bodies in a full 3D environment.
Since their size is fixed, and the elements are just 9,
they can be easily implemented without the need of
the blocks presented in the paper, but the use of these
blocks can make their implementation extremely fast
and straightforward.
The picture below shows the rotation matrix that
transforms body fixed coordinates into earth fixed
coordinates, given the orientation of the body with
respect to earth in roll, pitch and yaw angles [5].
psi
psi
Toc
Toc
1
Selector
m
theta
theta
ftp
Tcf
phi
phi
Tfb
Tfb
Z=X*Y
vmult1
Tcf
Z=X*Y
vmult2
2
Vb
Z=X*Y
1
vmult3
Ve
Figure 10. Body2Earth Rotation Matrix
The scheme above is clearly related to the formula:
1
0
0
0
1
0
e
Rb = e e e
ψS( 0 ) ϑS( 1 ) φS( 0 )
The Multiplication, Transposition, and Inversion
blocks, are currently being used by the first author in
a recursive frequency domain identification scheme
to be used for on-line flight parameter estimation and
self adaptive flight control.
Modularity issues
In our opinion, there’s also another fact, less evident
but very important, that makes a totally Simulinkbased implementation best than one which is
partially coded in Matlab or C.
Often, the persons who designs a Simulink scheme,
has barely taken a quick computer programming
course in his/her studies, so it’s pointless to expect
that he/she can write efficient code with good
modularity and reusability characteristics.
In these cases, a totally Simulink-based project
generally retains better modularity and reusability
characteristics than a partially coded project.
Since as project complexity grows modularity
becomes critical, this fact alone can, in our opinion,
become the most compelling reason to adopt a totally
Simulink-based implementation.
Where to find the blocks
The four blocks, along with all the examples in the
paper and instructions files for their use, can be
freely downloaded at the official mathworks ftp site:
ftp://ftp.mathworks.com/pub/contrib/v5/simulink/ma
trix_library
DSP Blockset
If you have the DSP toolbox, then you can use the
DSP matrix library, which is very efficient and
complete. The presented blocks are fully compatible
with the ones of the DSP matrix library, and of
course all the considerations above on the
opportunity of using the presented blocks, applies as
well to the blocks in the DSP matrix library.
Conclusions
In this paper, we briefly explained Simulink’s Sfunctions mechanism, which is used to represent and
simulate dynamical systems, in particular we have
shown the three common methods to build Sfunction, and underlined that in some cases neither
of them can be at the same time easy and suitable for
real time simulations.
Two Simulink blocks were proposed that solve this
problem by allowing varying matrices to be handled
directly in Simulink environment. Their aspect, their
use, and the most important part of their code were
explained, and some examples were given.
Finally some good reasons that render the use of
these blocks highly desirable were given.
References
[1] Matlab 5.3 and Simulink 3 manuals, The
Mathworks Inc.
[2] G. Campa, M. Sharma, A.J. Calise, M.
Innocenti, “Neural Network Augmentation of
Linear Controllers with Application to
Underwater Vehicles” American Control
Conference 2000, June 2-4, 2000 San Diego.
[3] F. Nardi, R.T. Rysdyk, A.J. Calise, “Neural
Network Based Adaptive Control of a Thrust
Vectored Ducted Fan,” AIAA Guidance,
Navigation and Control Conference, AIAA 993996, Portland, OR, 1999.
[4] F.L. Lewis, S. Jaganathan, A. Yesilderek,
Neural Network Control of Robot Manipulators
and Nonlinear Systems, Taylor & Francis,
London, 1999, pp. 150-170.
[5] Guidance and Control of Ocean Vehicles, Thor
J. Fossen John Wiley and Sons 1994.
6
American Institute of Aeronautics and Astronautics