Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.

php
ot100_granfm3.qxp
10/11/2007
9:29 AM
Page 1

Creating Simulations
with Simulink, Volume I
Numerical Computing
ot100_granfm3.qxp 10/11/2007 9:29 AM Page 3

Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Numerical Computing
with Simulink, Volume I
Creating Simulations

Richard J. Gran
Mathematical Analysis Company
Norfolk, Massachusetts

Society for Industrial and Applied Mathematics • Philadelphia

ot100_granfm3.qxp 10/11/2007 9:29 AM Page 4

Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Copyright © 2007 by the Society for Industrial and Applied Mathematics.

10 9 8 7 6 5 4 3 2 1
All rights reserved. Printed in the United States of America. No part of this book may be
reproduced, stored, or transmitted in any manner without the written permission of the
publisher. For information, write to the Society for Industrial and Applied Mathematics,
3600 Market Street, 6th floor, Philadelphia, PA 19104-2688 USA.
Trademarked names may be used in this book without the inclusion of a trademark symbol.
These names are used in an editorial context only; no infringement of trademark is intended.
Maple is a trademark of Maplesoft, Waterloo, Ontario, Canada.
MATLAB, Simulink, Real Time Workshop, SimHydraulics, Stateflow, and Handle Graphics are
registered trademarks of The MathWorks, Inc. For MATLAB product information, please
contact The MathWorks, Inc., 3 Apple Hill Drive, Natick, MA 01760-2098 USA, 508-647-
7000, Fax: 508-647-7101, info@mathworks.com, www.mathworks.com.
Figure 1.12 is used with permission of Marcus Orvando, president of Erwin Sattler Clocks
of America.
Figure 1.16 was taken by Michael Reeve on January 30, 2004, and its source is the
Wikipedia article “Foucault Pendulum.” Per Wikipedia, permission is granted to copy, dis-
tribute and/or modify this document under the terms of the GNU Free Documentation
License, Version 1.2 or any later version published by the Free Software Foundation; with no
Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. Subject to disclaimers.
Figure 2.7 is taken from the Wikipedia article “Thermostat.” Per Wikipedia, permission is
granted to copy, distribute and/or modify this document under the terms of the GNU Free
Documentation License, Version 1.2 or any later version published by the Free Software
Foundation; with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
Subject to disclaimers.

Library of Congress Cataloging-in-Publication Data

Gran, Richard J., 1940-
Numerical computing with Simulink / Richard J. Gran.
v. cm.
Includes bibliographical references and index.
Contents: v. 1. Creating simulations
ISBN 978-0-898716-37-5 (alk. paper)
1. Numerical analysis--Data processing. 2. SIMULINK. I. Title.

QA297.G676 2007
518.0285--dc22
2007061803

is a registered trademark.
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
ot100_granfm3.qxp
10/11/2007
9:29 AM
Page 5


To Dr. Richard A. Scheuing: A Friend and Mentor
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Contents

List of Figures xi

List of Tables xvii

Preface xix

1 Introduction to Simulink 1
1.1 Using a Picture to Write a Program . . . . . . . . . . . . . . . . . . . . . 1
1.2 Example 1: Galileo Drops Two Objects from the Leaning Tower of Pisa . 5
1.3 Example 2: Modeling a Pendulum and the Escapement of a Clock . . . . 17
1.3.1 History of Pendulum Clocks . . . . . . . . . . . . . . . . . . . . . 18
1.3.2 A Simulation Model for the Clock . . . . . . . . . . . . . . . . . . 20
1.4 Example 3: Complex Rotations—The Foucault Pendulum . . . . . . . . . 24
1.4.1 Forces from Rotations . . . . . . . . . . . . . . . . . . . . . . . . 25
1.4.2 Foucault Pendulum Dynamics . . . . . . . . . . . . . . . . . . . . 26
1.5 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2 Linear Differential Equations, Matrix Algebra, and Control Systems 33

2.1 Linear Differential Equations: Linear Algebra . . . . . . . . . . . . . . . 33
2.1.1 Solving a Differential Equation at Discrete Time Steps . . . . . . . 36
2.1.2 Linear Differential Equations in Simulink . . . . . . . . . . . . . . 38
2.2 Laplace Transforms for Linear Differential Equations . . . . . . . . . . . 40
2.3 Linear Feedback Control . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.3.1 What Is a Control System? . . . . . . . . . . . . . . . . . . . . . . 44
2.3.2 Control Systems and Linear Differential Equations . . . . . . . . . 48
2.4 Linearization and the Control of Linear Systems . . . . . . . . . . . . . . 49
2.4.1 Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
2.4.2 Eigenvalues and the Response of a Linear System . . . . . . . . . . 51
2.5 Poles and the Roots of the Characteristic Polynomial . . . . . . . . . . . . 55
2.5.1 Feedback of the Position of the Mass in the Spring-Mass Model . . 56
2.5.2 Feedback of the Velocity of the Mass in the Spring-Mass Model . . 57
2.5.3 Comparing Position and Rate Feedback . . . . . . . . . . . . . . . 61
2.5.4 The Structure of a Control System: Transfer Functions . . . . . . . 62

vii
 viii Contents
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

2.6 Transfer Functions: Bode Plots . . . . . . . . . . . . . . . . . . . . . . . 64

2.6.1 The Bode Plot for Continuous Time Systems . . . . . . . . . . . . 65
2.6.2 Calculating the Bode Plot for Continuous Time Systems . . . . . . 65
2.7 PD Control, PID Control, and Full State Feedback . . . . . . . . . . . . . 69
2.7.1 PD Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
2.7.2 PID Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.7.3 Full State Feedback . . . . . . . . . . . . . . . . . . . . . . . . . . 72
2.7.4 Getting Derivatives for PID Control or Full State Feedback . . . . 74
2.8 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

3 Nonlinear Differential Equations 81

3.1 The Lorenz Attractor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
3.1.1 Linear Operating Points: Why the Lorenz Attractor Is Chaotic . . . 83
3.2 Differential Equation Solvers in MATLAB and Simulink . . . . . . . . . 86
3.3 Tables, Interpolation, and Curve Fitting in Simulink . . . . . . . . . . . . 87
3.3.1 The Simple Lookup Table . . . . . . . . . . . . . . . . . . . . . . 88
3.3.2 Interpolation: Fitting a Polynomial to the Data and Using the Result
in Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
3.3.3 Using Real Data in the Model: From Workspace and File . . . . . . 92
3.4 Rotations in Three Dimensions: Euler Rotations, Axis-Angle
Representations, Direction Cosines, and the Quaternion . . . . . . . . . . 94
3.4.1 Euler angles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.4.2 Direction Cosines . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
3.4.3 Axis-Angle Rotations . . . . . . . . . . . . . . . . . . . . . . . . . 99
3.4.4 The Quaternion Representation . . . . . . . . . . . . . . . . . . . . 99
3.5 Modeling the Motion of a Satellite in Orbit . . . . . . . . . . . . . . . . . 105
3.5.1 Creating an Attitude Error When Using Direction Cosines . . . . . 107
3.5.2 Creating an Attitude Error Using Quaternion Representations . . . . 109
3.5.3 The Complete Spacecraft Model . . . . . . . . . . . . . . . . . . . 109
3.6 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

4 Digital Signal Processing in Simulink 115

4.1 Difference Equations, Fibonacci Numbers, and z-Transforms . . . . . . . 116
4.1.1 The z-Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
4.1.2 Fibonacci (Again) Using z-Transforms . . . . . . . . . . . . . . . . 120
4.2 Digital Sequences, Digital Filters, and Signal Processing . . . . . . . . . . 121
4.2.1 Digital Filters, Using z-Transforms, and Discrete Transfer Functions 121
4.2.2 Simulink Experiments: Filtering a Sinusoidal Signal and Aliasing . 123
4.2.3 The Simulink Digital Library . . . . . . . . . . . . . . . . . . . . . 128
4.3 Matrix Algebra and Discrete Systems . . . . . . . . . . . . . . . . . . . . 130
4.4 The Bode Plot for Discrete Time Systems . . . . . . . . . . . . . . . . . . 135
4.5 Digital Filter Design: Sampling Analog Signals, the Sampling Theorem,
and Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.5.1 Sampling and Reconstructing Analog Signals . . . . . . . . . . . . 137
 Contents ix
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

4.5.2 Analog Prototypes of Digital Filters: The Butterworth Filter . . . . 142

4.6 The Signal Processing Blockset . . . . . . . . . . . . . . . . . . . . . . . 145
4.6.1 Fundamentals of the Signal Processing Blockset: Analog Filters . . 146
4.6.2 Creating Digital Filters from Analog Filters . . . . . . . . . . . . . 148
4.6.3 Digital Signal Processing . . . . . . . . . . . . . . . . . . . . . . . 149
4.6.4 Implementing Digital Filters: Structures and Limited Precision . . . 153
4.6.5 Batch Filtering Operations, Buffers, and Frames . . . . . . . . . . . 160
4.7 The Phase-Locked Loop . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
4.8 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170
Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

5 Random Numbers, White Noise, and Stochastic Processes 173

5.1 Modeling with Random Variables in Simulink: Monte Carlo Simulations . 173
5.1.1 Monte Carlo Analysis and the Central Limit Theorem . . . . . . . . 174
5.1.2 Simulating a Rayleigh Distributed Random Variable . . . . . . . . 176
5.2 Stochastic Processes and White Noise . . . . . . . . . . . . . . . . . . . . 177
5.2.1 The Random Walk Process . . . . . . . . . . . . . . . . . . . . . . 178
5.2.2 Brownian Motion and White Noise . . . . . . . . . . . . . . . . . . 180
5.3 Simulating a System with White Noise Inputs Using the Weiner Process . 184
5.3.1 White Noise and a Spring-Mass-Damper System . . . . . . . . . . 184
5.3.2 Noisy Continuous and Discrete Time Systems:
The Covariance Matrix . . . . . . . . . . . . . . . . . . . . . . . . 186
5.3.3 Discrete Time Equivalent of a Continuous Stochastic Process . . . . 189
5.3.4 Modeling a Specified Power Spectral Density: 1/f Noise . . . . . . 194
5.4 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199

6 Modeling a Partial Differential Equation in Simulink 201

6.1 The Heat Equation: Partial Differential Equations in Simulink . . . . . . . 202
6.1.1 Finite Dimensional Models . . . . . . . . . . . . . . . . . . . . . . 202
6.1.2 An Electrical Analogy of the Heat Equation . . . . . . . . . . . . . 203
6.2 Converting the Finite Model into Equations for Simulation with Simulink 205
6.2.1 Using Kirchhoff’s Law to Get the Equations . . . . . . . . . . . . . 206
6.2.2 The State-Space Model . . . . . . . . . . . . . . . . . . . . . . . . 208
6.3 Partial Differential Equations for Vibration . . . . . . . . . . . . . . . . . 212
6.4 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213

7 Stateflow: A Tool for Creating and Coding State Diagrams, Complex Logic,
Event Driven Actions, and Finite State Machines 215
7.1 Properties of Stateflow: Building a Simple Model . . . . . . . . . . . . . 216
7.1.1 Stateflow Semantics . . . . . . . . . . . . . . . . . . . . . . . . . . 218
7.1.2 Making the Simple Stateflow Chart Do Something . . . . . . . . . 221
7.1.3 Following Stateflow’s Semantics Using the Debugger . . . . . . . . 223
7.2 Using Stateflow: A Controller for Home Heating . . . . . . . . . . . . . . 225
7.2.1 Creating a Model of the System and an Executable Specification . . 225
 x Contents
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

7.2.2 Stateflow’s Action Language Types . . . . . . . . . . . . . . . . . . 229

7.2.3 The Heating Controller Layout . . . . . . . . . . . . . . . . . . . . 230
7.2.4 Adding the User Actions, the Digital Clock, and the Stateflow Chart
to the Simulink Model of the Home Heating System . . . . . . . . . 231
7.2.5 Some Comments on Creating the GUI . . . . . . . . . . . . . . . . 239
7.3 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 240
Exercise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 240

8 Physical Modeling: SimPowerSystems and SimMechanics 241

8.1 SimPowerSystems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242
8.1.1 How the SimPowerSystems Blockset Works: Modeling a Nonlinear
Resistor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 246
8.1.2 Using the Nonlinear Resistor Block . . . . . . . . . . . . . . . . . 248
8.2 Modeling an Electric Train Moving on a Rail . . . . . . . . . . . . . . . . 251
8.3 SimMechanics: A Tool for Modeling Mechanical Linkages
and Mechanical Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 256
8.3.1 Modeling a Pendulum with SimMechanics . . . . . . . . . . . . . . 257
8.3.2 Modeling the Clock: Simulink and SimMechanics Together . . . . 260
8.4 More Complex Models in SimMechanics and SimPowerSystems . . . . . 262
8.5 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 265
Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 267

9 Putting Everything Together: Using Simulink in a System Design Process 269

9.1 Specifications Development and Capture . . . . . . . . . . . . . . . . . . 270
9.1.1 Modeling and Analysis: Converting the Specifications into an “Ex-
ecutable Specification” . . . . . . . . . . . . . . . . . . . . . . . . 271
9.2 Modeling the System to Incorporate the Specifications: Lunar Module
Rotation Using Time Optimal Control . . . . . . . . . . . . . . . . . . . 272
9.2.1 From Specification to Control Algorithm . . . . . . . . . . . . . . . 273
9.3 Design of System Components to Meet Specifications: Modify the Design
to Accommodate Computer Limitations . . . . . . . . . . . . . . . . . . . 276
9.3.1 Final Lunar Module Control System Executable Specification . . . 279
9.3.2 The Control System Logic: Using Stateflow . . . . . . . . . . . . . 283
9.4 Verification and Validation of the Design . . . . . . . . . . . . . . . . . . 285
9.5 The Final Step: Creating Embedded Code . . . . . . . . . . . . . . . . . 286
9.6 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 287

10 Conclusion: Thoughts about Broad-Based Knowledge 289

Bibliography 291

Index 295
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

List of Figures

1.1 Simulink library browser. . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Leaning Tower of Pisa. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Integrator block dialog. . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Simulink model for the Leaning Tower of Pisa experiment. . . . . . . . . 9
1.5 Scope block graph for the Leaning Tower simulation. . . . . . . . . . . . 10
1.6 Leaning Tower Simulink model with air drag added. . . . . . . . . . . . 13
1.7 Scope dialog showing how to change the number of axes in the Scope plot. 13
1.8 Leaning Tower simulation results when air drag is included. . . . . . . . 14
1.9 Adding the second object to the Leaning Tower simulation by vectorizing
the air drag. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.10 Simulating results. The Leaning Tower Simulink model with a heavy
and light object. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.11 Right clicking anywhere in the diagram opens a pull-down menu that
allows changes to the colors of the blocks, bold vector line styles, and
other model annotations. . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.12 Components of a pendulum clock. . . . . . . . . . . . . . . . . . . . . . 19
1.13 Simulink model of the clock. . . . . . . . . . . . . . . . . . . . . . . . . 21
1.14 Blocks in the subsystem “Escapement Model.” . . . . . . . . . . . . . . 22
1.15 Results of the clock simulation for times around 416 sec. Upper figure
is the pendulum angle, and the lower figure is the acceleration applied to
the pendulum by the escapement. . . . . . . . . . . . . . . . . . . . . . 24
1.16 Foucault pendulum at the Panthéon in Paris. . . . . . . . . . . . . . . . . 25
1.17 Axes used to model the Foucault pendulum. . . . . . . . . . . . . . . . . 27
1.18 The Simulink model of the Foucault pendulum. . . . . . . . . . . . . . . 29
1.19 Simulation results from the Foucault pendulum Simulink model using
the default solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
1.20 Simulation results for the Foucault pendulum using a tighter tolerance
for the solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.1 Using Simulink to compare discrete and continuous time state-space ver-
sions of the pendulum. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.2 Outputs as seen in the Scope blocks when simulating the model in Figure 2.1. 40
2.3 Linear pendulum model using Transfer Function and Zero-Pole-Gain
blocks from the Simulink Continuous Library. . . . . . . . . . . . . . . 43

xi
 xii List of Figures
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

2.4 Home heating system Simulink model (Thermo_NCS). This model de-
scribes the house temperature using a single first order differential equa-
tion and assumes that a bimetal thermostat controls the temperature. . . . 45
2.5 Subsystem “House” contains the single differential equation that models
the house. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
2.6 Results of simulating the home heating system model. The top figure
shows the outdoor and indoor temperatures, and the bottom figure shows
the cost of the heat over the 24-hour simulation. . . . . . . . . . . . . . . 47
2.7 A typical home-heating thermostat [50]. . . . . . . . . . . . . . . . . . . 47
2.8 Hysteresis curve for the thermostat in the Simulink model. . . . . . . . . 48
2.9 Model of a spring-mass-damper system using Simulink primitives. . . . 52
2.10 Model of a spring-mass-damper system using the state-space model and
Simulink’s automatic vectorization. . . . . . . . . . . . . . . . . . . . . 52
2.11 Changing the value of the damping ratio from 0.1 to 0.5 in the model of
Figure 2.10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
2.12 Generic components of a control system. . . . . . . . . . . . . . . . . . 55
2.13 A simple control that feeds back the position of the mass in a spring-
mass-damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
2.14 Feed back the velocity of the mass instead of the position to damp the
oscillations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
2.15 Time response of the mass (position at the top and velocity at the bottom)
for the velocity feedback controller. . . . . . . . . . . . . . . . . . . . . 58
2.16 Using the state-space model for the spring-mass-damper control system. . 59
2.17 Root locus for the spring-mass-damper system. Velocity gain varying
from 0 to 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
2.18 Root locus for the position feedback control. Gain changes only the
frequency of oscillation. . . . . . . . . . . . . . . . . . . . . . . . . . . 61
2.19 One method of determining the velocity of the mass is to differentiate the
position using a linear system that approximates the derivative. . . . . . 64
2.20 Simulation result from the Simulink model. . . . . . . . . . . . . . . . . 64
2.21 Control System Toolbox interface to Simulink. GUIs allow you to select
inputs and outputs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
2.22 LTI Viewer results from the LTI Viewer. . . . . . . . . . . . . . . . . . . 68
2.23 Using proportional-plus derivative (PD) control for the spring-mass-
damper. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
2.24 Response of the spring-mass-damper with the PD controller. There are
no oscillations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
2.25 Proportional, integral, and derivative (PID) control in the spring-mass-
damper system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
2.26 PID control of the spring-mass system. Response to a unit step has zero
error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
2.27 Getting speed from measurements of the mass position and force applied
to the mass. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
2.28 Simulink model to compare three methods for deriving the velocity of
the mass. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
2.29 Methods for deriving rate are selected using the “multiport switch.” . . . 77
 List of Figures xiii
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

3.1 Simulink model and results for the Lorenz attractor. . . . . . . . . . . . 83

3.2 Exploring the Lorenz attractor using the linearization tool from the control
toolbox. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
3.3 Eigenvalues of the linearized Lorenz attractor as a function of time. . . . 85
3.4 Pseudotabulated data for the external temperature in the home heating
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
3.5 Home heating system Simulink model with tabulated outside temperature
added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
3.6 Dialog to add the tabulated input. . . . . . . . . . . . . . . . . . . . . . 90
3.7 Real data must be in a MATLAB array for use in the “From Workspace”
block. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
3.8 Simulink model of the home heating system using measured outdoor
temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.9 From workspace dialog block that uses the measured outdoor temperature
array. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.10 Using real data from measurements in the Simulink model. . . . . . . . . 94
3.11 Single axis rotation in three dimensions. . . . . . . . . . . . . . . . . . . 95
3.12 Using quaternions to compute the attitude from three simultaneous body
axis rotations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
3.13 Simulink block to convert quaternions into a direction cosine. . . . . . . 103
3.14 Combining the quaternion block with the body rotational acceleration. . . 104
3.15 Electrical circuit of a motor with the mechanical equations. . . . . . . . . 106
3.16 Controlling the rotation of a spacecraft using reaction wheels. . . . . . . 110
3.17 Matrix concatenation blocks create the matrix Q used to achieve a desired
value for q. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

4.1 Simulating the Fibonacci sequence in Simulink. . . . . . . . . . . . . . . 117

4.2 Fibonacci sequence graph generated by the Simulink model. . . . . . . . 118
4.3 Simulink model for computing the golden ratio from the Fibonacci
sequence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
4.4 Generating the data needed to compute a digital filter transfer function. . 123
4.5 Adding callbacks to a Simulink model uses Simulink’s model properties
dialog, an option under the Simulink window’s file menu. . . . . . . . . 125
4.6 Sampling a sine wave at two different rates, illustrating the effect of
aliasing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
4.7 Changing the numerator and denominator in the Digital Filter block
changes the filter icon. . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
4.8 Using the Digital Filter block to simulate the Fibonacci sequence. . . . . 130
4.9 State-space models for discrete time simulations in Simulink. . . . . . . 131
4.10 A numerical experiment in Simulink: Does fk+1 fk−1 − fk2 = ±1? . . . . 134
4.11 13 iterations of the Fibonacci sequence show that fk+1 fk−1 − fk2 = ±1
(so far). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
4.12 The Bode plot for the discrete filter model developed using the Control
System Toolbox interface with Simulink. . . . . . . . . . . . . . . . . . 137
4.13 Illustrating the steps in the proof of the sampling theorem. . . . . . . . . 139
 xiv List of Figures
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

4.14 Simulation illustrating the sampling theorem. The input signal is 15

sinusoidal signals with frequencies less than 500 Hz. . . . . . . . . . . . 140
4.15 Specification of a unity gain low pass filter requires four pieces of infor-
mation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
4.16 Butterworth filter for D/A conversion using the results of the M-file “but-
terworthncs.” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
4.17 Using the Signal Processing Blockset to design an analog filter. . . . . . 147
4.18 Moving average FIR filter simulation. . . . . . . . . . . . . . . . . . . . 150
4.19 Pole-zero plot and filter properties for the moving average filter. . . . . . 152
4.20 The Simulink model for the band-pass filter with no computational limi-
tations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
4.21 The band-pass filter design created by the Signal Processing Blockset
digital filter block. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
4.22 Band-pass filter simulation results (with and without signal in the pass
band). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
4.23 Fixed-point implementation of a band-pass filter using the Signal Pro-
cessing Blockset. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
4.24 Running the Simulink model to determine the minima, maxima, and
scaling for all of the fixed-point calculations in the Band-pass filter. . . . 161
4.25 Illustrating the use of buffers. Reconstruction of a sampled signal using
the FFT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
4.26 Using the FFT and the sampling theorem to interpolate a faster sampled
version of a sampled analog signal. . . . . . . . . . . . . . . . . . . . . 163
4.27 The analog signal (sampled at 1 kHz, top) and the reconstructed signal
(sampled at 8 kHz, bottom). . . . . . . . . . . . . . . . . . . . . . . . . 165
4.28 Results of interpolating a sampled signal using FFTs and the sampling
theorem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
4.29 Phase-locked loop (PLL). A nonlinear feedback control for tracking fre-
quency and phase of a sinusoidal signal. . . . . . . . . . . . . . . . . . . 167
4.30 Voltage controlled oscillator subsystem in the PLL model. . . . . . . . . 168
4.31 Details of the integrator modulo 1 subsystem in Figure 4.30. . . . . . . . 168
4.32 Phase-locked loop simulation results. . . . . . . . . . . . . . . . . . . . 169

5.1 Monte Carlo simulation demonstrating the central limit theorem. . . . . . 175
5.2 Result of the Monte Carlo simulation of the central limit theorem. . . . . 175
5.3 A Simulink model that generates Rayleigh random variables. . . . . . . . 177
5.4 Simulink model that generates nine samples of a random walk. . . . . . . 179
5.5 Nine samples of a random walk process. . . . . . . . . . . . . . . . . . . 179
5.6 Simulink model for the spring-mass-damper system. . . . . . . . . . . . 185
5.7 Motions of 10 masses with white noise force excitations. . . . . . . . . . 185
5.8 White noise block in the Simulink library (masked subsystem). . . . . . 186
5.9 Continuous linear system covariance matrix calculation in Simulink. The
result is the covariance matrix that can be used for simulating with a
covariance equivalent discrete system. . . . . . . . . . . . . . . . . . . . 191
5.10 Noise response continuous time simulation and equivalent discrete sys-
tems at three different sample times. . . . . . . . . . . . . . . . . . . . . 192
 List of Figures xv
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

5.11 Simulation results from the four simulations with white noise inputs. . . 193
5.12 Simulation of the fractal noise process that has a spectrum proportional
to 1/f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196
5.13 Power spectral density estimators in the signal processing blockset used
to compute the sampled PSD of the 1/f noise process. . . . . . . . . . . 198

6.1 An electrical model of the thermodynamics of a house and its heating

system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204
6.2 Simulink model for the two-room house dynamics. . . . . . . . . . . . . 210
6.3 Simulation results when the heat is off. . . . . . . . . . . . . . . . . . . 211
6.4 Simulation results when the heat is on continuously. . . . . . . . . . . . 211
6.5 Using PID control to maintain constant room temperatures (70 deg F). . . 211

7.1 First Stateflow diagram (in the chart) uses a manual switch to create the
event. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 218
7.2 A Stateflow diagram that contains most of the stateflow semantics. . . . . 219
7.3 The Simulink model for the simple timer example. . . . . . . . . . . . . 222
7.4 Modified Stateflow chart provides the timer outputs “Start,” “On,” and
“Trip.” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223
7.5 Outputs resulting from running the simple timer example. . . . . . . . . 224
7.6 PID controller for the home-heating example developed in Chapter 6. . . 227
7.7 Response of the home heating system with the external temperature from
Chapter 3 as input. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228
7.8 A tentative design for the new home heating controller we want to
develop. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228
7.9 The complete simulation of the home heating controller. . . . . . . . . . 231
7.10 Graphical user interface (GUI) that implements the prototype design for
the controller. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232
7.11 Interactions between the GUI and Stateflow use changes in the gain values
in three gain blocks in the subsystem called “User_Selections.” . . . . . 233
7.12 The first of the two parallel states in the state chart for the heating
controller. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235
7.13 The second of the two parallel states in the home heating controller state
chart. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 236
7.14 The final version of the PID control for the home heating controller. . . . 239

8.1 A simple electrical circuit model using Simulink and SimPowerSystems. 242
8.2 The electrical circuit time response. Current in the resistor and voltage
across the capacitor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243
8.3 The SimPowerSystems library in the Simulink browser and the elements
sublibrary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 244
8.4 A Simulink and SimPowerSystems subsystem model for a nonlinear or
time varying resistor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 247
8.5 Masking the subsystem with an icon that portrays what the variable re-
sistor does makes it easy to identify in an electrical circuit. . . . . . . . . 248
8.6 Circuit diagram of a simulation with the nonlinear resistor block. . . . . 249
 xvi List of Figures
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

8.7 Using the signal builder block in Simulink to define the values for the
time varying resistance. . . . . . . . . . . . . . . . . . . . . . . . . . . 250
8.8 Using SimPowerSystems: Results from the simple RL network
simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251
8.9 Simulink and SimPowerSystems model of a train and its dynamics. . . . 252
8.10 Complete model of the train moving along a track. . . . . . . . . . . . . 254
8.11 The Simulink blocks that compute the rail resistance as a function of the
train location. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255
8.12 Simulation results. Rail resistances, currents in the rails, and current to
the train. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255
8.13 Modeling a pendulum using SimMechanics. . . . . . . . . . . . . . . . . 258
8.14 Coordinate system for the pendulum. . . . . . . . . . . . . . . . . . . . 258
8.15 Simulating the pendulum motion using SimMechanics and its viewer. . . 260
8.16 Adding Simulink blocks to a SimMechanics model to simulate the clock
as in Chapter 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261
8.17 Time history of the pendulum motion from the clock simulation. . . . . . 261
8.18 SimMechanics model Sim_Mechanics_Vibration. . . . . . . . . 263
8.19 String model (the subsystem in the model above). It consists of 20 iden-
tical “spring mass” revolute and prismatic joint elements. . . . . . . . . . 263
8.20 String simulation. Plucked 15 cm from the left, at the center, and 15 cm
from the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263
8.21 Using SimPowerSystems to build a model for a four-room house. . . . . 264

9.1 The start of the specification capture process. Gather existing designs
and simulations and create empty subsystems for pieces that need to be
developed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 271
9.2 The physics of the phase plane logic and the actions that need to be
developed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 276
9.3 Complete simulation of the lunar module digital autopilot. . . . . . . . . 279
9.4 Simulating a counter that ticks at 625 microsec for the lunar module
Simulink model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 281
9.5 The blocks in the jet-on-time counter subsystem in Figure 9.4. . . . . . . 281
9.6 Lunar module graph of the switch curves and the phase plane motion for
the yaw axis. (The graph uses a MATLAB function block in the simulation.)282
9.7 Stateflow logic that determines the location of the state in the phase
plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284
9.8 Simulink blocks that check different model attributes for model
verification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 286
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

List of Tables

3.1 Computation time for the Lorenz model with different solvers and toler-
ances. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.1 Using low pass filters. Simulating various filters with various sample
times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.2 Comparison of the error in converting a digital signal to an analog signal
using different low pass filters. . . . . . . . . . . . . . . . . . . . . . . . 148
4.3 Specification for the band-pass filter and the actual values achieved in the
design. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154

5.1 Simulink computation times for the four white noise simulations. . . . . 193

xvii
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Preface

Simulation has gradually become the preferred method for designing complex systems. The
aerospace industry provided the impetus for this design approach because most aerospace
systems are at best difficult or at worst impossible to test. I started my career working on the
lunar module flight control system. Testing this system in space was extremely expensive
and very difficult. In the entire development cycle, there was, in fact, only one unmanned test
in earth orbit. All of the designs, by necessity, required simulations. Over time, this approach
has become pervasive in many industries, particularly the automotive industry. Even in
product developments as mundane (seemingly) as disk drives for computers, simulation has
become the dominant method of design.
A digital simulation of the environment that a system “lives” in is crucial in this
modern simulation based design process. The accurate representation of this environment
allows the designer to see how well the system being designed performs. It allows the
designer to verify that the design meets all of the performance specifications. Moreover,
as the simulation develops, the designer always can see how the performance requirements
pass through to the various subsystems in the design. This refinement of the specification
allows the model to, in essence, become the specification. The simulation environment can
also allow the designer of embedded digital signal processing and control algorithms to test
the code directly as it is developed to ensure that it satisfies the specification. The last,
and newest, step in this process is the ability to use the simulation model automatically to
generate the computer code for computers in the system. This can eliminate most of the
“hand coding” that is required by the current practice. (In most companies this process
consists of handing a stack of papers with the written specification to a team of computer
programmers for coding.) This approach eliminates two potential sources of error: the
errors introduced by the conversion of the design into a written specification and the errors
introduced by the manual coding.
Simulink® is a remarkable tool that fills the niche for a robust, accurate, and easily used
simulation tool. It provides a visual way of constructing a simulation of a complex system
that evolves in time (the designer can optimize the system by optimizing the mathematical
model’s response). This “model-based design” approach provides accurate simulation as
various numerical constants in the model are changed. It allows the user to tune these
constants so that the simulation accurately portrays the real world. It then allows sharing
of the model among many different users so they can begin to design various components
or parts of the system. At the same time, it provides a means for understanding design
problems whose solutions will make the system operate in a desired way. Finally, the
resulting design’s embedded software is rapidly convertible into embedded code that is

xix
 xx Preface
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

directly usable for testing in a stand-alone computer. All of this can be accomplished using
one easily understood graphical environment. A major thesis of this book is that Simulink
can form the basis for an integrated design environment. As your use of the tool evolves,
the way this can happen should become more apparent.
Because simulation is such an important adjunct to numerical analysis, the book can
be a supplement for a course on numerical computing. It fits very well with an earlier
SIAM book, Numerical Computing with MATLAB by Cleve Moler [29] (which we refer to
so frequently that its title is abbreviated NCM in the text), and it has been written with the
same eclectic mix of practice and mathematics that makes Cleve’s book so much fun to read.
I hope that this work favorably compares with—and comes close to matching the quality
of—Cleve’s work. As is the case with Numerical Computing with MATLAB, this book is
based on a collection of 88 model files that are in the library Numerical Computing with
Simulink (which we abbreviate as NCS), which can be downloaded from the Mathematical
Analysis Co. web page (www.math-analysis.com). This library of SIMULINK models and
MATLAB® M-files forms an essential part of the book and, through extensions that the
reader is encouraged to peruse, an essential part of the exercises.
The book will take you on a tour of the Simulink environment, showing you how to
develop a system model and then how to execute the design steps to make the model into
a functioning design laboratory. Along the way, you will be introduced to the mathematics
of systems, including difference equations and z-transforms, ordinary differential equations
(both linear and nonlinear), Laplace transforms, numerical methods for solving differential
equations, and methods for simulating complex systems from several different disciplines.
The mathematics of simulation is not complete without a discussion of random vari-
ables and random processes for doing Monte Carlo simulations. Toward this end, we
introduce and develop the techniques for modeling random processes with predetermined
statistical properties. The mathematics for this type of simulation begins with “white noise,”
which is a very difficult entity to pin down. We introduce the concept of Brownian motion
and show the connection between this process and the fictitious, but useful, white noise
process. The simulation of the Brownian process in Simulink is developed; it is then used
to form more (statistically) complex processes. We review and show how to simulate and
analyze random processes, and we formulate and use the power spectral density of a process.
We introduce other tools, in addition to Simulink, from The MathWorks. Each of
the tools is an expansion into a different “domain of knowledge.” The first tool is the
Signal Processing Blockset that extends Simulink into the domain of signal processing
(both analog and digital). The second tool, Stateflow® , expands Simulink to include state
charts and signal flow for modeling event driven systems (i.e., systems where the actions
start at times that are asynchronous with the computer’s clock). This tool naturally extends
Simulink into the traditional realm of finite state machines.
The third tool, SimPowerSystems, extends Simulink into the realm of physical mod-
eling and in particular into the realm of electrical circuits including power systems, motor
drives, power generation equipment including power electronics, and three-phase power
transmission lines.
The last tool, SimMechanics, is also a physical modeling tool; it develops models of
mechanical systems such as robots, linkages, etc. These tools can all work together in the
common Simulink environment, but they all have their own method for displaying a picture
of the underlying system.
 Preface xxi
Downloaded 09/26/12 to 128.104.1.219. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Both MATLAB and Simulink are available in student versions. All of the examples
in this book use the tools in the student version supplemented with the Control System
Toolbox, the Signal Processing blockset, SimMechanics, and SimPowerSystems, which
can be purchased separately. If you are not a student but wish to learn how to use Simulink
using this book, the MathWorks can provide a demonstration copy of MATLAB, Simulink,
and the add-ons; contact them at http://www.mathworks.com.
This book not only introduces Simulink and discusses many useful tricks to help
develop models but it also shows how it is possible for Simulink to be the basis of a
complete design process. I hope the material is not too theoretical and that you find the
models we will work with are fun, instructive, and useful in the future.
As is always the case, innumerable people have contributed to the genesis and devel-
opment of this book. I thank Cleve Moler for sharing an early manuscript of his wonderful
book and The MathWorks for many years of interesting and stimulating work. In particu-
lar, I wish to thank the MathWorks staff with whom I have had many long and interesting
discussions: Rob Aberg, Paul Barnard, Jason Ghidella, Ned Gully, Loren Shure, and Mark
Ullman, and The MathWorks’s CEO, Jack Little. These interesting and productive discus-
sions helped me create the various examples and models described in this book.
Words are never enough to thank ones family. My wife, and soulmate, Jean Frova
Gran, provides support, joy, affection, and love. My daughter Kathleen keeps my thoughts
focused on the future. Her hard work and dedication to the achievement of her goals, despite
many small-business frustrations, is laudable.
My extended family includes two stepgrandchildren, Ian and Ila, who provide me
with insight into the joy brought by the process of learning about this unbelievable world
we live in. My stepdaughters Elizabeth and Juliet have a sense of humor and a joie de vivre
that are inspirational.
Finally, I need to say “thank you” to my sister Cora and her wonderful husband,
Arthur, who have always been loving and supportive in more ways than can be stated; Cora
is, as the self-proclaimed “matriarch” of the Grans, much more to me than a mere sister.

Dr. Richard J. Gran

Mathematical Analysis Co.
www.math-analysis.com
Norfolk, Massachusetts, August 2006