Using Simulink in Engine Model PDF

S IMULINK -S TATEFLOW T ECHNICAL E XAMPLES
Using Simulink and Stateflow

® TM
in Automotive Applications
Abstract This book includes nine examples that represent typical design tasks of an automotive engineer. It
shows how The MathWorks modeling and simulation tools, Simulink® and Stateflow, facilitate
TM
the design of automotive control systems. Each example explains the principles of the physical sit-
uation, and presents the equations that represent the system. The examples show how to proceed
from the physical equations to the Simulink block diagram. Once the Simulink model has been
completed, we run the simulation, analyze the results, and draw conclusions from the study.
TABLE OF C ONTENTS
Introduction ....................................................................................................................... 4
System Models in Simulink................................................................................................ 7

I. Engine Model.............................................................................................................. 8
II. Anti-Lock Braking System........................................................................................ 18
III. Clutch Engagement Model....................................................................................... 23
IV. Suspension System.................................................................................................... 31
V. Hydraulic Systems .................................................................................................... 35
System Models in Simulink with Stateflow Enhancements .......................................... 49

VI. Fault-Tolerant Fuel Control System ........................................................................ 50
VII. Automatic Transmission Control............................................................................. 61
VIII. Electrohydraulic Servo Control ............................................................................... 71
IX. Modeling Stick-Slip Friction .................................................................................... 84
USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 3

INTRODUCTION
Summary Automotive engineers have found simulation to be a vital tool in the timely and cost-effective
development of advanced control systems. As a design tool, Simulink has become the standard for
excellence through its flexible and accurate modeling and simulation capabilities. As a result of its open
architecture, Simulink allows engineers to create custom block libraries so they can leverage each other’s
work. By sharing a common set of tools and libraries, engineers can work together effectively within
individual work groups and throughout the entire engineering department.
In addition to the efficiencies achieved by Simulink, the design process can also benefit from Stateflow, an
interactive design tool that enables the modeling and simulation of complex reactive systems. Tightly
integrated with Simulink, Stateflow allows engineers to design embedded control systems by giving them
an efficient graphical technique to incorporate complex control and supervisory logic within their
Simulink models.
This booklet describes nine automotive design examples that illustrate the strengths of Simulink and
Stateflow in accelerating and facilitating the design process.
Examples The examples cited in this booklet consist of application design tasks typically encountered in
Description the automotive industry. We present a variety of detailed models including the underlying
equations, block diagrams, and simulation results. The material may serve as a starting point for the
new Simulink user or as a reference for the more experienced user. In the models, we propose approaches
for model development, present solutions to challenging problems, and illustrate some of the most
common design uses of Simulink and Stateflow today.
The applications and models described in this booklet include the following examples using Simulink
alone:
I. Engine Model
engine.mdl — open-loop simulation
enginewc.mdl — closed-loop simulation
II. Anti-Lock Braking System
absbrake.mdl
III. Clutch Engagement Model
clutch.mdl
IV. Suspension System
suspn.mdl
V. Hydraulic Systems
hydcyl.mdl — Pump and actuator assembly
hydcyl4.mdl — Four-cylinder model
hydrod.mdl — Two-cylinder model with load constraints
4 S IMULINK -S TATEFLOW T ECHNICAL EXAMPLES

The following applications and models use Simulink enhanced with Stateflow:
VI. Fault-Tolerant Fuel Control System
fuelsys.mdl
VII. Automatic Transmission Control

sf_car.mdl
VIII. Electrohydraulic Servo Control

sf_electrohydraulic.mdl
IX. Modeling Stick-Slip Friction

sf_stickslip.mdl
Simulink The models used in this book are available via ftp at
Model Files ftp://ftp.mathworks.com/pub/product-info/examples/autobook.zip. This zip file contains the set
of MDL, MAT, and M-files containing Simulink models that users can explore and build upon. The
included files require MATLAB® 5.1, Simulink 2.1, and Stateflow 1.0. Models for these applications can be
opened in Simulink by typing the name of the model at the MATLAB command prompt. MATLAB,
Simulink, and Stateflow are not included with this booklet. To obtain a copy of MATLAB, Simulink, and
Stateflow, or for a diskette containing the model files, please contact your representative at The
MathWorks.
Acknowledgments The engine model is based on published findings by Crossley and Cook (1991)(1). We’d like to thank
Ken Butts and Jeff Cook of the Ford Motor Company for permission to include this model and for
subsequent help in building the model in Simulink.
The clutch and hydraulic cylinder models are based on equations provided by General Motors. We’d like
to thank Eric Gassenfeit of General Motors for permission to include these models.
The vehicle suspension model was written by David MacClay of Cambridge Control Ltd.
The simple three-state engine model and the set of icons that are relevant for automotive modeling were
provided by Modular Systems. A far more detailed engine model may be purchased directly from Modular
Systems.
Contact The MathWorks technical personnel specializing in automotive solutions can be reached via e-mail
Information at the following addresses:
Stan Quinn squinn@mathworks.com
Andy Grace agrace@mathworks.com
Paul Barnard pbarnard@mathworks.com
Larry Michaels lmichaels@mathworks.com
Bill Aldrich baldrich@mathworks.com

Or contact any of our international distributors and resellers directly. See the back page for additional
contact information.
Both Modular Systems and Cambridge Control Ltd. offer consulting services in automotive modeling.
They can be reached as follows:
Attention: Robert W. Weeks

Modular Systems
714 Sheridan Road
Evanston, IL 60202-2502 USA
Tel: 708-869-2023
E-mail: bobweeks@ix.netcom.com
Attention: Sham Ahmed

Cambridge Control Ltd.
Newton House
Cambridge Business Park
Cowley Road
Cambridge, DB4 4WZ UK
011/44-1223-423-2
E-mail: Sham@camcontrol.co.uk

System Models in Simulink

I. ENGINE MODEL
Summary This example presents a model of a four-cylinder spark ignition engine and demonstrates Simulink’s
capabilities to model an internal combustion engine from the throttle to the crankshaft output. We used
well-defined physical principles supplemented, where appropriate, with empirical relationships that
describe the system’s dynamic behavior without introducing unnecessary complexity.
Overview This example describes the concepts and details surrounding the creation of engine models with emphasis
on important Simulink modeling techniques. The basic model uses the enhanced capabilities of
Simulink 2 to capture time-based events with high fidelity. Within this simulation, a triggered
subsystem models the transfer of the air-fuel mixture from the intake manifold to the cylinders via
discrete valve events. This takes place concurrently with the continuous-time processes of intake flow,
torque generation and acceleration. A second model adds an additional triggered subsystem that provides
closed-loop engine speed control via a throttle actuator.
These models can be used as standalone engine simulations. Or, they can be used within a larger system
model, such as an integrated vehicle and powertrain simulation, in the development of a traction control
system.
Model Description
This model, based on published results by Crossley and Cook (1991), describes the simulation of a four-
cylinder spark ignition internal combustion engine. The Crossley and Cook work also shows how a
simulation based on this model was validated against dynamometer test data.
The ensuing sections (listed below) analyze the key elements of the engine model that were identified by
Crossley and Cook:
• Throttle
• Intake manifold
• Mass flow rate
• Compression stroke
• Torque generation and acceleration
Note: Additional components can be added to the model to provide greater accuracy in simulation and to
more closely replicate the behavior of the system.
Analysis THROTTLE
and Physics The first element of the simulation is the throttle body. Here, the control input is the angle of the throttle
plate. The rate at which the model introduces air into the intake manifold can be expressed as the product
of two functions—one, an empirical function of the throttle plate angle only; and the other, a function of
the atmospheric and manifold pressures. In cases of lower manifold pressure (greater vacuum), the flow
rate through the throttle body is sonic and is only a function of the throttle angle. This model accounts for

this low pressure behavior with a switching condition in the compressibility equations shown in
Equation 1.1.
m˙ ai = f (θ )g (Pm )
= mass flow rate into manifold (g/s) where,
f (θ ) = 2.821 − 0.05231θ + 0.10299θ 2 − 0.00063θ 3
θ = throttle angle (deg) Equation 1.1
 P
1, Pm ≤ amb
 2
 2 2 P
 Pm Pamb − Pm , amb
≤ Pm ≤ Pamb
g (Pm ) =  Pamb 2
 2 2
− Pm Pamb − Pamb , Pamb ≤ Pm ≤ 2Pamb
 m P
−1, Pm ≥ 2Pamb

Pm = manifold pressure (bar)
Pamb = ambient (atmospheric) pressure (bar)
Intake Manifold
The simulation models the intake manifold as a differential equation for the manifold pressure. The
difference in the incoming and outgoing mass flow rates represents the net rate of change of air mass with
respect to time. This quantity, according to the ideal gas law, is proportional to the time derivative of the
manifold pressure. Note that, unlike the model of Crossley and Cook, 1991(1) (see also references 3
through 5), this model doesn’t incorporate exhaust gas recirculation (EGR), although this can easily be
added.
(m˙ ai − m˙ ao )
RT
P˙ m = Equation 1.2
Vm
where,
R = specific gas constant

T = temerature (˚K)
Vm = manifold volume (m 3 )
m˙ ao = mass flow rate of air out of the manifold (g/s)
P˙ m = rate of change of manifold pressure (bar/s)
Intake Mass Flow Rate

The mass flow rate of air that the model pumps into the cylinders from the manifold is described in
Equation 1.3 by an empirically derived equation. This mass rate is a function of the manifold pressure
and the engine speed.
2
m˙ ao = −0.366 + 0.08979 NPm − 0.0337NPm + 0.0001N 2P m Equation 1.3

where,
N = engine speed (rad/s)

Pm = manifold pressure (bar)
To determine the total air charge pumped into the cylinders, the simulation integrates the mass flow rate
from the intake manifold and samples it at the end of each intake stroke event. This determines the total
air mass that is present in each cylinder after the intake stroke and before compression.
Compression Stroke
In an inline four-cylinder four-stroke engine, 180° of crankshaft revolution separate the ignition of each
successive cylinder. This results in each cylinder firing on every other crank revolution. In this model,
the intake, compression, combustion, and exhaust strokes occur simultaneously (at any given time, one
cylinder is in each phase). To account for compression, the combustion of each intake charge is delayed
by 180° of crank rotation from the end of the intake stroke.
Torque Generation and Acceleration

The final element of the simulation describes the torque developed by the engine. An empirical
relationship dependent upon the mass of the air charge, the air/fuel mixture ratio, the spark advance, and
the engine speed is used for the torque computation.
Torque eng = −181.3 + 379.36 m a + 21.91( A / F ) − 0.85( A / F ) 2 + 0.26σ − 0.0028σ 2

Equation 1.4
+0.027N − 0.000107N 2 + 0.00048Nσ + 2.55σm a − 0.05σ 2 m a
where,
ma = mass of air in cylinder for combustion (g)

A / F = air to fuel ratio
σ = spark advance (degrees before top - dead - center
Torqueeng = torque produced by the engine (Nm)
The engine torque less the net load torque results in acceleration.
JN˙ = Torqueeng − Torqueload Equation 1.5
where,
J = Engine rotational moment of inertia (kg-m2)

Ṅ = Engine acceleration (rad/s2)

Modeling — We incorporated the model elements described above into an engine model using Simulink. The
The Open-Loop following sections describe the decisions we made for this implementation and the key Simulink elements
Simulation used. This section shows how to implement a complex nonlinear engine model easily and quickly in the
Simulink environment. We developed this model in conjunction with Ken Butts, Ford Motor Company (2).
Figure 1.1 shows the top level of the Simulink model. Note that, in general, the major blocks correspond
to the high-level list of functions given in the model description in the preceding summary. Taking
advantage of Simulink’s hierarchical modeling capabilities, most of the blocks in Figure 1.1 are made up
of smaller blocks. The following paragraphs describe these smaller blocks.
Engine Timing Model in Simulink 2

choose Start from A Demonstration of Triggered Subsystems
the Simulation
menu to run
edge180 N 1
crank speed
valve timing (rad/sec)
mass(k) Air Charge
Torque Teng
mass(k+1)
Throttle Ang. 30/pi
N N
trigger
throttle Mass Airflow Rate
(degrees) Tload
1 rad/s Engine
Engine Speed, N Combustion to
Compression Speed
s rpm
Vehicle (rpm)
Throttle & Manifold Dynamics
Load
Intake
Drag Torque throttle deg (purple)
load torque Nm (yellow)
Mux
Figure 1.1: The top level of the Simulink engine model
Throttle/Manifold
Simulink models for the throttle and intake manifold subsystems are shown in Figure 1.2. The throttle
valve behaves in a nonlinear manner and is modeled as a subsystem with three inputs. Simulink
implements the individual equations, given in Equation 1.1 as function blocks. These provide a
convenient way to describe a nonlinear equation of several variables. A Switch block determines whether
the flow is sonic by comparing the pressure ratio to its switch threshold, which is set at one half (Equation
1.1). In the sonic regime, the flow rate is a function of the throttle position only. The direction of flow is
from the higher to lower pressure, as determined by the Sign block. With this in mind, the Min block
ensures that the pressure ratio is always unity or less.

The intake manifold is modeled by the differential equation as described in Equation 1.2 to compute the
manifold pressure. A Simulink function block also computes the mass flow rate into the cylinder, a
function of manifold pressure and engine speed (Equation 1.3).
Throttle Manifold Dynamics

1 Throttle Angle, theta (deg)
Throttle Ang. Limit to Positive

Manifold Pressure, Pm (bar) Throttle Flow, mdot (g/s)
1.0 Atmospheric Pressure, Pa (bar)
Atmospheric Throttle
Pressure, Pa
(bar)
mdot Input (g/s) mdot to Cylinder (g/s) 1

Mass Airflow Rate
2 N (rad/sec) Manifold Pressure, Pm (bar)
Engine Speed, N
Intake Manifold
Throttle Flow vs. Valve Angle and Pressure
Throttle Angle,
theta (deg) f(theta)
1 2.821 - 0.05231*u + 0.10299*u*u - 0.00063*u*u*u
2
Manifold Pressure, g(pratio)
Pm (bar)
min pratio 2*sqrt(u - u*u)
1
3
1.0 Throttle
Atmospheric Pressure, Flow, mdot
Pa (bar) Sonic Flow (g/s)
flow direction
Intake Manifold Vacuum

2
Manifold Pressure,
Pm (bar)
1
1
mdot Input 0.41328
(g/s) s
RT/Vm Mu -0.366 + 0.08979*u[1]*u[2] - 0.0337*u[2]*u[1]*u[1] + 0.0001 *u[1]*u[2]*u[2] 1
p0 = 0.543 bar
Pumping mdot to
2 Cylinder
(g/s)
N (rad/sec)
Figure 1.2: The Throttle and Intake Manifold Subsystems

Intake and Compression
An integrator accumulates the cylinder mass air flow in the Intake block. The Valve Timing block issues
pulses that correspond to specific rotational positions in order to manage the intake and compression
timing. Valve events occur each cam rotation, or every 180° of crankshaft rotation. Each event triggers a
single execution of the Compression subsystem. The output of the trigger block within the Compression
subsystem then feeds back to reset the Intake integrator. In this way, although both triggers conceptually
occur at the same instant in time, the integrator output is processed by the Compression block immediately
prior to being reset. Functionally, the Compression subsystem uses a Unit Delay block to insert 180° (one
event period) of delay between the intake and combustion of each air charge.
Consider a complete four-stroke cycle for one cylinder. During the intake stroke, the Intake block
integrates the mass flow rate from the manifold. After 180° of crank rotation, the intake valve closes and
the Unit Delay block in the Compression subsystem samples the integrator state. This value, the
accumulated mass charge, is available at the output of the Compression subsystem 180° later for use in
combustion. During the combustion stroke, the crank accelerates due to the generated torque. The final
180°, the exhaust stroke, ends with a reset of the Intake integrator, prepared for the next complete 720°
cycle of this particular cylinder.
For four cylinders, we could use four Intake blocks, four Compression subsystems, etc., but each would be
idle 75% of the time. We’ve made the implementation more efficient by performing the tasks of all four
cylinders with one set of blocks. This is possible because, at the level of detail we’ve modeled, each function
applies to only one cylinder at a time.
Combustion
Engine torque is a function of four variables. The model uses a Mux block to combine these variables into
a vector that provides input to the Torque Gen block. Here, a function block computes the engine torque,
as described empirically in Equation 1.4. The torque which loads the engine, computed by step functions
in the Drag Torque block, is subtracted in the Vehicle Dynamics subsystem. The difference divided by the
inertia yields the acceleration, which is integrated to arrive at the engine crankshaft speed.
Results We saved the Simulink model in the file engine.mdl which can be opened by typing engine at the
MATLAB prompt. Select Start from the Simulation menu to begin the simulation. Simulink scope
windows show the engine speed, the throttle commands which drive the simulation, and the load torque
which disturbs it. Try adjusting the throttle to compensate for the load torque.
Figure 1.3 shows the simulated engine speed for the default inputs:
 8.97, t < 5
Throttle(deg) = 
11.93, t ≥ 5

25, t≤ 2

Load ( Nm ) =  20, 2 < <8
t
 25, t≥ 8

Engine Simulation
3500
3000
2500
engine speed (RPM)
2000
1500
1000
500
0
0 1 2 3 4 5 6 7 8 9 10
time (sec)
Figure 1.3: The simulated engine speed
Note the behavior as the throttle angle and load torque change.
Modeling — SPEED C ONTROL
The Closed-Loop The following enhanced model demonstrates the flexibility and extensibility of Simulink models. In the
Simulation enhanced model, the objective of the controller is to regulate engine speed with a fast throttle actuator,
such that changes in load torque have minimal effect. This is easily accomplished in Simulink by adding
a discrete-time PI controller to the engine model as shown in Figure 1.4.
The model is stored in the file enginewc.mdl, which can be opened by typing enginewc at the MATLAB
command prompt. This represents the same engine model described previously but with closed-loop
control.

choose Start from
the Simulation
Closed Loop Engine Speed Control
menu to run
edge180 N 1
crank speed
valve timing (rad/sec)
mass(k) Air Charge

Desired rpm
Torque Teng
mass(k+1)
speed Throttle Setting Throttle Ang.
N N 30/pi
set trigger
N
point Mass Airflow Rate
Tload
1 rad/s Engine
Engine Speed, N Combustion to
Controller Compression Speed
s rpm
Vehicle (rpm)
Throttle & Manifold Dynamics
Load
Intake
drag torque
throttle deg (purple)
load torque Nm (yellow)
Mux
Figure 1.4: A discrete-time PI controller is added to the engine model

to regulate speed
We chose a control law which uses proportional plus integral (PI) control. The integrator is needed to
adjust the steady-state throttle as the operating point changes, and the proportional term compensates for
phase lag introduced by the integrator.
∫
θ = K p ( N set − N ) + K I ( N set − N ) dt ,
N set = speed set point Equation 1.6
K p = proportional gain
K I = integral gain
A discrete-time controller, suitable for microprocessor implementation, is employed. The integral term in
Equation 1.6 must thus be realized with a discrete-time approximation.
As is typical in the industry, the controller execution is synchronized with the engine’s crankshaft
rotation. The controller is embedded in a triggered subsystem that is triggered by the valve timing signal
described above. The detailed construction of the Controller subsystem is illustrated in Figure 1.5. Of note
is the use of the Discrete-Time Integrator block with its sample time parameter set (internally) at -1. This
indicates that the block should inherit its sample time, in this case executing each time the subsystem is
triggered. The key component that makes this a triggered subsystem is the Trigger block shown at the
bottom of Figure 1.5. Any subsystem can be converted to a triggered subsystem by dragging a copy of this
block into the subsystem diagram from the Simulink Connections library.

Triggered PI Controller
1 pi/30
Desired Kp
rpm rpm
to Proportional Gain 1
rad/s Ki T Throttle Setting
limit
z -1 output
Integral Gain
2 N
Discrete Time
0 Integrator
integrator input
enable integration
controller output
prevent windup
Figure 1.5: Speed Controller Subsystem
Results Typical simulation results are shown in Figure 1.6. The speed set point steps from 2000 to 3000 RPM
at t = 5 sec. The torque disturbances are identical to those used in the previous example. Note the quick
transient response, with zero steady-state error.
Several alternative controller tunings are shown. These can be adjusted by the user at the MATLAB
command line. This allows the engineer to understand the relative effects of parameter variations.
Closed–Loop Engine Speed Control

3200
3000
Kp = 0.05, Ki = 0.1
Kp = 0.033, Ki = 0.064
Kp = 0.061, Ki = 0.072
2800
engine speed (RPM)
2600
2400
2200
2000
1800
0 1 2 3 4 5 6 7 8 9 10
time (Sec.)
Figure 1.6: Typical simulation results

Conclusions The ability to model nonlinear, complex systems, such as the engine model described here, is one of
Simulink’s key features. The power of the simulation is evident in the presentation of the models above.
Simulink retains model fidelity, including precisely timed cylinder intake events, which is critical in
creating a model of this type. The two different models, the basic engine and complete speed control
system, demonstrate the flexibility of Simulink. In particular, the Simulink modeling approaches allow
rapid prototyping of an interrupt-driven engine speed controller.
References
1. P.R. Crossley and J.A. Cook, IEE International Conference ‘Control 91’, Conference Publication 332,
vol. 2, pp. 921-925, 25-28 March, 1991, Edinburgh, U.K.
2. The Simulink Model. Developed by Ken Butts, Ford Motor Company. Modified by Paul Barnard, Ted
Liefeld and Stan Quinn, The MathWorks, Inc., 1994-7.
3. J. J. Moskwa and J. K. Hedrick, "Automotive Engine Modeling for Real Time Control Application,"
Proc.1987 ACC, pp. 341-346.
4. B. K. Powell and J. A. Cook, "Nonlinear Low Frequency Phenomenological Engine Modeling and
Analysis," Proc. 1987 ACC, pp. 332-340.
5. R. W. Weeks and J. J. Moskwa, "Automotive Engine Modeling for Real-Time Control Using
Matlab/Simulink," 1995 SAE Intl. Cong. paper 950417.

There are two inputs to the Vehicle Suspension model shown in Figure 4.2. The first input is the road
height. A step input here corresponds to the vehicle driving over a road surface with a step change in
height. The second input is a horizontal force acting through the center of the wheels that results from
braking or acceleration maneuvers. Since the longitudinal body motion is not modeled, this input
appears only as a moment about the pitch axis.
The spring/damper subsystem that models the front and rear suspensions is shown in Figure 4.3. The
block is used to model Equation 4.1 through 4.3. The equations are implemented directly in the Simulink
diagram through the straightforward use of Gain and Summation blocks. The differences between front
and rear are accounted for as follows. Because the subsystem is a masked block, a different data set (L, K
and C) can be entered for each instance. Furthermore, L is thought of as the Cartesian coordinate x, being
negative or positive with respect to the origin, or center of gravity. Thus, Kf , Cf and -Lf are used for the
front and Kr, Cr and Lr for the rear.
L 1
L
2*K pitch
MomentArm3 Torque
MomentArm1
stiffness
2
L
Vertical
Force
MomentArm2 2*C
1 em
THETA
THETAdot damping
Z Fz
Zdot
Two DOF Spring/Damper Model
Figure 4.3: The Spring/Damper suspension subsystem
Results To run this model, first set up the required parameters in the MATLAB workspace. Run the following
M-file by typing suspdat, or from the MATLAB command line, enter the data by typing:
Lf = 0.9; % front hub displacement from body CG

Lr = 1.2; % rear hub displacement from body CG
Mb = 1200; % body mass in kg
Iyy = 2100; % body moment of inertia about y-axis in kgm^2
kf = 28000; % front suspension stiffness in N/m
kr = 21000; % rear suspension stiffness in N/m
cf = 2500; % front suspension damping in N/(m/s)
cr = 2000; % rear suspension damping in N/(m/s)
To run the simulation, select Start from the Simulink Simulation menu or type the following at the
MATLAB command line:
[t,x] = sim('suspn2',10); % run a time response

Figure 4.4 shows the plotted output results. You can automate setting the parameters, running the
simulation, and plotting these graphs by typing suspgrph at the MATLAB command line prompt.
-3 Vehicle Suspension Model Simulation

x 10
5
THETAdot
dθ/dt
0
-5
0.1
dz/dt
Zdot
-0.1
7000
Ff
6500 reaction force at front wheels
6000 -3
x 10
15
10
5 road height
h
0
-5
100
moment due to vehicle accel/decel
50
Y
0
0 1 2 3 4 5 6 7 8 9 10
time in seconds
Figure 4.4: A summary of the suspension simulation output results
Conclusions The Vehicle Suspension model allows you to simulate the effects of changing the suspension damping
and stiffness, thereby investigating the tradeoff between comfort and performance. In general, a racing
car has very stiff springs with a high damping factor, whereas a passenger vehicle designed for comfort
has softer springs and a more oscillatory response.

V. HYDRAULIC S YSTEMS
Summary This example considers several hydraulic systems. The general concepts apply to suspension, brake,
steering, and transmission systems. We model three variations of systems employing pumps, valves, and
cylinder/piston actuators. The first features a single hydraulic cylinder which we develop, simulate and
save as a library block. In the next model, we use four instances of this block, as in an active suspension
system. In the final model, we model the interconnection of two hydraulic actuators, held together by a
rigid rod which supports a large mass.
In some cases we treat relatively small volumes of fluid as incompressible. This results in a system of
differential-algebraic equations (DAEs). Simulink solvers are well-suited to handle this type of problem
efficiently. The masking and library reference capabilities add extra power and flexibility. The creation of
custom blocks enables the implementation of important subsystems with user-defined parameter sets.
The Simulink library keeps a master version of these blocks so that models using a master block
automatically incorporate any revisions and refinements made to it.
Analysis and Figure 5.1 shows a schematic diagram of the basic model. The model directs the pump flow Q to
Physics supply pressure p1 from which laminar flow q1ex leaks to exhaust. The control valve for the
piston/cylinder assembly is modeled as turbulent flow through a variable-area orifice. Its flow q12 leads to
intermediate pressure p2 which undergoes a subsequent pressure drop in the line connecting it to the
actuator cylinder. The cylinder pressure p3 moves the piston against a spring load, resulting in position x.
control valve
Q p1 p2
q23
q1ex A q12 C1
p3
pump
v3
C2
cylinder x
Figure 5.1: Schematic diagram of the basic hydraulic system
At the pump output, the flow is split between leakage and flow to the control valve. The leakage, q1ex is
modeled as laminar flow.

Q = q12 + q1ex
q1ex = C 2 p1
p1 = ( Q − q12 ) /C 2
Q = pump flow
Equation 5.1
q12 = control valve flow
q1ex = leakage
C 2 = flow coefficient
p1 = pump pressure
We modeled turbulent flow through the control valve with the orifice equation. The sign and absolute
value functions accommodate flow in either direction.
2
q12 = C d A sgn( p1 − p 2 ) p1 − p 2
ρ
C d = orifice discharge coefficient
A = orifice area Equation 5.2
p 2 = pressure downstream of control valve
ρ = fluid density
The fluid within the cylinder pressurizes due to this flow, q12 = q23, less the compliance of the piston
motion. We also modeled fluid compressibility in this case.
β
p˙ 3 =
v3
(q 12 − xA
˙ c )
p 3 = piston pressure
β = fluid bulk modulus
v 3 = fluid volume at p 3 Equation 5.3
=V 30 + Ac x
Ac = cylinder
cylinder cross–sectional
cross – sectional area
area
V 30 = fluid volume at x = 0
We neglected the piston and spring masses due to the large hydraulic forces. Force balance at the piston
gives:
x = p3 Ac / K
Equation 5.4
K = spring rate

We complete the system of equations by differentiating this relationship and incorporating the pressure
drop between p2 and p3. The latter models laminar flow in the line from the valve to the actuator.
x˙ = p˙ 3 A c /K
q23 = q12
= C 1( p 2 − p 3 ) Equation 5.5
p 2 = p 3 + q12 /C 1
C 1 = laminar flow coefficient
Modeling Figure 5.2 shows the basic model, stored in the file hydcyl.mdl. Simulation inputs are the pump
flow and the control valve orifice area. The model is organized as two subsystems — the pump and the
actuator assembly.
Hydraulic Cylinder Model
Mux
Qout p1 p
p1 pressures
p1 (yellow)
x
p2 (purple)
A p3 (blue)
pump qin piston position
valve/cylinder/piston/spring assembly
control valve
orifice area
Double click to run the Double click to see

Simulation for 0.1 seconds a 4 cylinder model.
Figure 5.2: The basic pump/valve/actuator model
Pump
The pump model computes the supply pressure as a function of the pump flow and the load (output) flow
(Figure 5.3). A From Workspace block provides the pump flow data, Qpump. This is specified by a matrix
with column vectors of time points and the corresponding flow rates [T, Q]. The model subtracts the
output flow, using the difference, the leakage flow, to determine the pressure p1, as indicated above in
Equation 5.1. Since Qout = q12 is a direct function of p1 (via the control valve), this forms an algebraic
loop. An estimate of the initial value, p10, enables a more efficient solution.
We mask the subsystem in Simulink to facilitate ready access to the parameters by the user. The
parameters to be specified are T, Q, p10 and the leakage flow coefficient, C2. For easy identification, we then
assigned the masked block the icon shown in Figure 5.2, and saved it in the Simulink library hydlib.mdl.

IC
[T,Q]
[p10]
1/C2 1
Qpump 1
p1
Qout leakage
Figure 5.3: Hydraulic Pump Subsystem
Actuator Assembly
valve/cylinder/piston/spring assembly
p2, p3
1/C1
p2
xdotAc laminar flow
Ac^2/K pressure drop Mux 1
p
p2
A beta
2 Area q12
1 p3 Ac/K 2
1 p1 s
V30 x
p1 control valve piston
force/spring
flow pressure
cylinder
volume
3
Ac
qin
Figure 5.4: Hydraulic Actuator Subsystem
In Figure 5.4, a system of differential-algebraic equations models the cylinder pressurization with the
pressure p3, which appears as a derivative in Equation 5.3 and is used as the state (integrator). If we
neglect mass, the spring force and piston position are direct multiples of p3 and the velocity is a direct
multiple of ṗ 3 . This latter relationship forms an algebraic loop around the bulk modulus Gain block,
Beta . The intermediate pressure p2 is the sum of p3 and the pressure drop due to the flow from the valve to
the cylinder (Equation 5.5). This relationship also imposes an algebraic constraint through the control
valve and the 1/C1 gain.
The control valve subsystem computes the orifice (equation 5.2), with the upstream and downstream
pressures as inputs, as well as the variable orifice area. A lower level subsystem computes the “signed
square root,” y = sgn(u) u . Three nonlinear functions are used, two of which are discontinuous. In
combination, however, y is a continuous function of u.

Results B ASELINE M ODEL
We simulated the model with the following data.
Cd = 0.61
ρ = 800 kg/m3
C 1 = 2e- 8 m3 /sec/Pa
C 2 = 3e- 9 m3 /sec/Pa
β = 7e8 Pa
Ac = 0.001 m3
K = 5e4 N/m
V30 = 2.5e- 5 m3
We specified the pump flow as:
sec. m3/sec
 0 0.005
 
 0.04 0.005
 0.04 0 
[T, Q] =  
 0.05 0 
 0.05 0.005
 
 0.1 0.005
The system thus initially steps to a pump flow of 0.005 m3/sec = 300 l/min, abruptly steps to zero at t =
0.04, then resumes its initial flow rate at t = 0.05.
The control valve starts with zero orifice area and ramps to 1e-4 m2 during the 0.1 second simulation
time. With the valve closed, all of the pump flow goes to leakage so the initial pump pressure jumps to
p 10 = Q/C2 = 1667 KPa.
As the valve opens, pressures p2 and p3 build up while p1 dips in response to the load increase as shown in
Figure 5.5. When the pump flow cuts off, the spring and piston act like an accumulator and p3, though
decreasing, is continuous. The flow reverses direction, so p2, though relatively close to p 3, falls abruptly.
At the pump itself, all of the backflow goes to leakage and p 1 drops radically. This behavior reverses as the
flow is restored.

5 Baseline Hydraulic Pressure Response
x 10
18
16
p1
14 p1
12
system pressures (Pa)

p2
10
p3
8
p2
p
3
p
2
6
4 p3 p
1
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time (sec)
Figure 5.5: Pressures in baseline simulation
The piston position is directly proportional to p3 , where the hydraulic and spring forces balance as shown
in Figure 5.6. Discontinuities in the velocity at 0.04 and 0.05 seconds indicate negligible mass. The model
reaches a steady state when all of the pump flow again goes to leakage, now due to zero pressure drop
across the control valve. That is, p3 = p2 = p1 = p10.
Baseline Hydraulic Cylider

0.035
0.03
0.025
piston position (m)
0.02
0.015
0.01
0.005
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time (sec)
Figure 5.6: Baseline hydraulic cy linder p iston p osition

We have now tested the pump and actuator blocks and determined that they perform according to design.
We created a Simulink library and copied these blocks into it. After saving the library (hydlib.mdl), we
replaced the blocks in the original model with library copies. The model file, hydcyl.mdl, now contains
references to the library blocks rather than all of the details of the subsystems. This demonstrates how we
build master libraries of important system components. Other designers can now employ identical copies
of these blocks in other systems. Whenever improvements are made to the library blocks, Simulink
automatically incorporates the changes into each individual model.
Four Cylinder Model

We now construct a new model with a single pump and four actuators (Figure 5.7). The same pump
pressure p1 drives each cylinder assembly and the sum of their flows loads the pump. Although each of
the four control valves could be controlled independently, as in an active suspension system, in this case
all four receive the same commands, a linear ramp in orifice area from zero to 0.002 m2.
Four Cylinder Model
p
p1
x
A
qin
pump
valve/cylinder/piston assembly 1
p positions
p1
x1 (yellow)
x x2 (purple)
A x3 (blue)
qin x4 (red)
control valve/cylinder/piston assembly 2
Mux
valve
command
p
p1
x
A
qin
p
p1
x
A
qin
supply
pressure
p1 load
flow
Double click for a model

Double click to run the with two hydraulic cylinders
Simulation for 0.1 seconds interconnected by a rigid rod
Figure 5.7: A single pump d riving four actuators
The pump flow begins at 0.005 m3/sec again for this system, then drops to half that value at t = 0.05 sec.
The parameters C1, C2, Cd, ρ and V30 are identical to the previous model. However, by assigning

individual values for K, Ac and, in one case, β, each of the four cylinders exhibit distinctive transient
responses. In relationship to the parameter values used above, the model characterizes the four actuators
according to Table 5.1.
parameter actuator 1 actuator 2 actuator 3 actuator 4
spring rate K K/4 4K K

piston area Ac Ac/4 4Ac Ac
bulk modulus β β β β/1000
Table 5.1: Parameter c omparison for individual actuators
The ratio of area/spring rate remains constant, so each case should have the same steady-state output.
The dominant time constant for each subsystem is proportional to Ac 2/K, so we can expect case two to be
somewhat faster than the case one, and case three somewhat slower. In case four, the effective bulk
modulus of the fluid is significantly lower, as would be the case with entrained air. We thus expect this
softer case to respond more sluggishly than case one. The simulation results support these predictions.
Piston Positions in Four Cylinder Example

0.02
0.018
0.016
0.014
piston position (m)
0.012
0.01
x2
0.008
x1
0.006
0.004 x4
0.002
x3
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time (sec)
Figure 5.8: Actuator Positions for Four-Cylinder Model
The initial jolt of flow at t = 0 responds like a pressure impulse, as seen by the four actuators (in Figure
5.8). The pump pressure p1 , which is initially high, drops rapidly as all four loads combine to make a
high flow demand. During the initial transient (about four milliseconds), distinctive responses identify
the individual dynamic characteristics of each unit.

As predicted by the differences in parameter values, actuator two responds much faster than the baseline,
actuator one. The third and fourth devices are much slower because they require more oil to move the
same distance. In case three, the piston displaces more volume due to its larger cross-sectional area. In
case four, although the displaced volume is the same as in case one, the device requires more oil because it
is subsequently compressed.
The distinctions in behavior are blurred, however, as the pump pressure falls to the level within the
cylinders (in Figure 5.9). The individual responses blend into an overall system response which
maintains the flow balance between the components. At t = 0.05 seconds, the pump flow drops to a level
that is close to equilibrium and the actuator flows are nearly zero. The individual steady-state piston
positions are equal, as predicted by design.
-4 Flow Transients in Four Cylinder Example

x 10
20
q
3
15
q4
actuator flow rate (m /sec)
10
3
5
q
1
q
2
-5
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time (sec)
Figure 5:9: Individual flow rates in four cylinder model
Two Cylinder Model with Load Constraints

In the final model (Figure 5.10), a rigid rod which supports a large mass interconnects two hydraulic
actuators. The model eliminates the springs as it applies the piston forces directly to the load. These forces
balance the gravitational force and result in both linear and rotational displacement.

Two Cylinder Model with Connecting Rod
Mux
z (purple), za (blue)
& zb (yellow)
theta
(rad)
zb
za
zbdot
theta
zadot
Mechanical
Load
Fext
Fb
Fa
- 9.81*M
F
p
qin
qin
valve/cylider/piston valve/cylider/piston
force assembly 1 force assembly 2
xdot
xdot
p1
p1
A
A
x
x
orifice B orifice A
supply pressure
load
flow pump
Double click to run the

Simulation for 0.1 seconds
Figure 5.10: Two hydraulic c ylinders with c onnecting rod
The load subsystem shown in Figure 5.11 solves the equations of motion, which we compute directly with
standard Simulink blocks. We assume the rotation angle is small.
Mz˙˙ = F b + F a + F ext
z = displacement at center
M = total mass
F a , F b = piston forces
F ext = external force at center Equation 5.6
L L
Iθ˙˙ = Fb − Fa
2 2
θ = angular displacement, clockwise
I = moment of inertia
L = rod length

The positions and velocities of the individual pistons follow directly from the geometry.
L
za = z − θ
2
L
zb =z + θ
2
L ˙
z˙ a = z˙ − θ Equation 5.7
2
L
z˙ b = z˙ + θ˙
2
z a ,z b = piston displacements
Connecting Rod
2
Fext
1 zdot 1
1 1/M 3
s s
Fb z
zdot zbdot 2 position

zbdot translation
thetadot zadot 6 z zb 1
zadot zb
velocity
translation theta za 5
za
1 1
L/2 1/I 4
s thetadot s
theta
3
Fa
Figure 5.11: Mechanical load subsystem
The parameters used in the simulation are identical to the first model, except:
L = 1.5m
M = 2500 kg
I = 100 kg- m2
Q pump = 0.005 m3 /sec (constant)
C2 = 3e- 10 m3 /sec/Pa
F ext = -9.81M N
Although pump flow is constant in this case, the model controls the valves independently according to the
following schedule in Figure 5.12:

-5 System Control Valve Inputs
x 10
1.2
B A
1
control valve orifice area (m2)

0.8
0.6
0.4
0.2
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time (secs)
Figure 5.12: Complementary a ctuator c ontrol action
Figure 5.13 and Figure 5.14 show the simulation outputs of rod displacement and angle, respectively. The
response of z is typical of a type-one (integrating) system. The relative positions and the angular
movement of the rod illustrate the response of the individual actuators to their out-of-phase control
signals.

-3 Piston and Load Positions
x 10
12
10
linear piston displacement (m)

z
a
6
z
z
b
4
-2
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time (sec)
Figure 5.13: Linear piston and load motion
-4 Two Cylinder Model Rotational Response

x 10
4
0
load angular displacement (rad)
-2
-4
-6
-8
-10
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time (sec)
Figure 5.14: Rotational displacement of load

Conclusions Simulink provides a productive environment for simulating hydraulic systems, offering enhancements
that provide enormous productivity in modeling and flexibility in numerical methods. The use of masked
subsystems and model libraries facilitates structured modeling with automatic component updates. That
is, as users modify library elements, the models that use the elements automatically incorporate the new
versions. Simulink can use differential-algebraic equations (DAEs) to model some fluid elements as
incompressible and others as compliant, allowing efficient solutions for complex systems of
interdependent circuits.
Models such as these will ultimately be used as part of overall plant or vehicle systems. The hierarchical
nature of Simulink allows independently developed hydraulic actuators to be placed, as appropriate, in
larger system models, for example, adding controls in the form of sensors or valves. In cases such as
these, tools from the MATLAB Control System Toolbox can analyze and tune the overall closed-loop
system. The MATLAB/Simulink environment can thus support the entire design, analysis, and modeling
cycle.

System Models in Simulink with Stateflow Enhancements

VI. FAULT -T OLERANT FUEL CONTROL S YSTEM
Summary The following example illustrates how to combine Stateflow with Simulink to efficiently model hybrid
systems. This type of modeling is particularly useful for systems that have numerous possible
operational modes based on discrete events. Traditional signal flow is handled in Simulink while
changes in control configuration are implemented in Stateflow.
The model described below represents a fuel control system for a gasoline engine. The system is highly
robust in that individual sensor failures are detected and the control system is dynamically reconfigured
for uninterrupted operation.
Analysis and Similar to the engine model described earlier in this document, physical and empirical relationships
Physics form the basis for the throttle and intake manifold dynamics of this model. The mass flow rate of air
pumped from the intake manifold, divided by the fuel rate which is injected at the valves, gives the air-fuel
ratio. The ideal, or stoichiometric mixture ratio provides a good compromise between power, fuel
economy, and emissions. A target ratio of 14.6 is assumed in this system.
Typically, a sensor determines the amount of residual oxygen present in the exhaust gas (EGO). This
gives a good indication of the mixture ratio and provides a feedback measurement for closed-loop control.
If the sensor indicates a high oxygen level, the control law increases the fuel rate. When the sensor detects
a fuel-rich mixture, corresponding to a very low level of residual oxygen, the controller decreases the fuel
rate.
Modeling Figure 6.1 shows the top level of the Simulink model (fuelsys.mdl). The controller uses signals from
the system’s sensors to determine the fuel rate which gives a stoichiometric mixture. The fuel rate
combines with the actual air flow in the engine gas dynamics model to determine the resulting mixture
ratio as sensed at the exhaust.
The user can selectively disable each of the four sensors (throttle angle, speed, EGO and manifold absolute
pressure [MAP]), to simulate failures. Simulink accomplishes this with Manual Switch blocks. Double-
click on the block itself to change the position of the switch. Similarly, the user can induce the failure
condition of a high engine speed by toggling the switch on the far left. A Repeating Table block provides the
throttle angle input and periodically repeats the sequence of data specified in the mask.

Fault Tolerant Fuel Control System
Choose Start from

the Simulation throttle
menu to run command
the model. engine speed o2_out
throttle
0 throttle
sensor throttle angle MAP
Nominal engine
Speed speed
engine speed
300
700
0 speed
High fuel rate fuel air/fuel ratio
Speed sensor
(rad./Sec.)
EGO
engine
gas
Use this switch
12 dynamics
to force the EGO
engine to sensor
overspeed
MAP
Metered Fuel
To toggle 0 MAP
a switch, sensor fuel rate
double click controller
on its icon. air/fuel
Use these switches to
simulate any combination mixture ratio
of sensor failures
Figure 6.1: Simulink fuelsys model
The controller uses the sensor input and feedback signals to adjust the fuel rate to give a stoichiometric
ratio (Figure 6.2). The model uses four subsystems to implement this strategy: control logic, sensor
correction, airflow calculation, and fuel calculation. Under normal operation, the model estimates the
airflow rate and multiplies the estimate by the reciprocal of the desired ratio to give the fuel rate. Feedback
from the oxygen sensor provides a closed-loop adjustment of the rate estimation in order to maintain the
ideal mixture ratio.
1
Sensor correction and
throtle Fault Redundancy
Airflow calculation
Mu Sensors
Corrected sens_in
Failures
Fuel Calculation
2
engine
Failures airflow airflow
speed
control logic fuel rate 1
fuel
mode
throt throtState mode rate
3 speedState
speed Mu
EGO
o2State
Ego
pressState
fuel_mode
4 press
MAP
fuel rate controller
Figure 6.2: Fuel rate c ontroller subsystem

Control Logic
A single Stateflow chart, consisting of a set of six parallel states, implements the control logic in its
entirety. The four parallel states shown at the top of Figure 6.3 correspond to the four individual sensors.
The remaining two parallel states at the bottom consider the status of the four sensors simultaneously
and determine the overall system operating mode. The model synchronously calls the entire Stateflow
diagram at a regular sample time interval of 0.01 sec. This permits the conditions for transitions to the
correct mode to be tested on a timely basis.
fuel rate
controller/control logic
Oxygen_Sensor_Mode O2_normal
Pressure_Sensor_Mode
entry: o2State = 0
[t > o2_t_thresh]
[press > max_press | press < min_press]
O2_warmup /Sens_Failure_Counter.INC
entry: o2State = 1 [Ego > max_ego]/
Sens_Failure_Counter.INC press_norm press_fail
entry: pressState = 0 entry: pressState = 1
[Ego < max_ego] / O2_fail

Sens_Failure_Counter.DEC entry: o2State = 1 [press > min_press & press < max_press] /
Sens_Failure_Counter.DEC
Throttle_Sensor_Mode
Speed_Sensor_Mode
[throt> max_throt | throt < min_throt]/
Sens_Failure_Counter.INC [speed==0 & press < zero_thresh]/
Sens_Failure_Counter.INC
throt_norm throt_fail
entry: throtState=0 entry: throtState = 1
speed_norm speed_fail
entry: speedState = 0 entry: speedState = 1
[throt > min_throt & throt < max_throt] [speed > 0] /

/ Sens_Failure_Counter.DEC Sens_Failure_Counter.DEC
Sens_Failure_Counter
MultiFail
INC
INC INC INC
FL0 FL1 FL2 FL4

FL3
DEC DEC DEC DEC
Fueling_Mode
Fuel_Disabled
H [ speed > max_speed ] entry: fuel_mode = DISABLED
Running
Overspeed
Low_Emmisions Rich_Mixture
entry: fuel_mode = LOW H entry: fuel_mode = RICH
[in(FL1)]
Single_Failure [!in(MultiFail)] [in(speed_norm) & ...
Normal
speed < (max_speed - hys)]
[in(FL0)]
enter(MultiFail) [in(MultiFail)]
[in(FL1)]
Warmup Shutdown
[in(O2_normal)]
exit(MultiFail)
Printed 14 Nov 1997 12:53:27
Figure 6.3: Control logic S tateflow diagram

When execution begins, all of the states start in their “normal” mode with the exception of the oxygen
sensor. The O2_warmup state is entered initially until time has exceeded the o2_t_thresh. The system
detects throttle and pressure sensor failures when their measured values fall outside their nominal
ranges. A manifold vacuum in the absence of a speed signal indicates a speed sensor failure. The oxygen
sensor also has a nominal range for failure conditions but, because zero is both the minimum signal
level and the bottom of the range, failure can be detected only when it exceeds the upper limit.
Regardless of which sensor fails, the model always generates the directed event broadcast
Sens_Failure_Counter.INC. In this way the triggering of the universal sensor failure logic is
independent of the sensor. The model also uses a corresponding sensor recovery event,
Sens_Failure_Counter.DEC. The Sens_Failure_Counter state keeps track of the number of failed
sensors. The counter increments on each Sens_Failure_Counter.INC event and decrements on each
Sens_Failure_Counter.DEC event. The model uses a superstate, MultiFail, to group all cases where
more than one sensor has failed.
The bottom parallel state represents the fueling mode of the engine. If a single sensor fails, operation
continues but the air/fuel mixture is richer to allow smoother running at the cost of higher emissions. If
more than one sensor has failed, the engine shuts down as a safety measure, since the air/fuel ratio
cannot be controlled reliably.
During the oxygen sensor warm-up, the model maintains the mixture at normal levels. If this is
unsatisfactory, the user can change the design by moving the warm-up state to within the Rich_Mixture
superstate. If a sensor failure occurs during the warm-up period, the Single_Failure state is entered
after the warm-up time elapses. Otherwise, the Normal state is activated at this time.
We easily added a protective overspeed feature by creating a new state in the Fuel_Disabled superstate.
Through the use of history junctions, we assured that the chart returns to the appropriate state when the
model exits the overspeed state. As the safety requirements for the engine become better specified, we can
add additional shutdown states to the Fuel_Disabled superstate.
Sensor Correction
The Fault Correction block determines which sensors to use and which to estimate. Figure 6.4 shows the
block diagram for this subsystem. The failures input is a vector of logic signals that trigger the
application of estimates to each particular sensor. When a component of the signal is nonzero, it enables
the appropriate estimation subsystem and causes the switch relating to that signal to send the estimate as
the output. Since the estimation routines are within enabled subsystems, they do not introduce any
computational overhead when they are not needed.

Throttle Estimate
Sensors throtle
throttle Speed Estimate
Sensors we
speed
1
m Corrected
EGO MAP Estimate
Sensors map
MAP
throttle
sensor
failure
1
Sensors
2 m speed sensor failure
Failures
pressure sensor failure
6.4: Sensor c orrection and fault redundancy
The sensors input to the Correction block is the vector of raw sensor values. When there are no faults, the
input simply passes on as the output signal. When a fault exists, the appropriate estimation block uses
this signal to recover the missing component. Figure 6.5 shows an estimation example of the algorithm
for the manifold pressure sensor.
MAP Estimation
speed
1
Sensors
1
map
throttle
Pressure Estimate (2D)
Enable
Figure 6.5: Manifold a bsolute p ressure reconstruction

Airflow Calculation
The Airflow Calculation block (Figure 6.6) is the location for the central control laws. The block estimates
the intake air flow to determine the fuel rate which gives the appropriate air/fuel ratio. Closed-loop
control adjusts the estimation according to the residual oxygen feedback in order to maintain the mixture
ratio precisely. Even when a sensor failure mandates open-loop operation, the most recent closed-loop
adjustment is retained to best meet the control objectives.
Intake Airflow Estimation and Closed Loop Correction

Throttle transient
correction
0.01z-0.01
throttle angle z - 0.8
1
est.
engine speed, N air
flow
1 em
sens_in
Pumping Constant
manifold pressure, Pm
Feedforward Control
0.5
Ramp
Rate (Ki)
EGO, residual
exhaust oxygen
e1
<= e0 T e2
2
0.5 z -1
feedback
Oxygen Sensor Integrator correction
Switching Threshold hold
0 integrator
2
O2 fail
Failures (warmup) Feedback Control
NOR
3 enable integration
~=
mode LOW not normal operation
Figure 6.6: Airflow Estimation and Correction
The engine’s intake air flow can be formulated as the product of the engine speed, the manifold pressure
and a time-varying scale factor.

N P
q= V cd η m
4π RT
= C pump( N , P m ) NP m ,
= intake mass flow
N = engine speed (rad/sec)
V cd = engine cylinder displacement volume Equation 6.1
η = volumetric efficiency
P m = manifold pressure
R = specific gas constant
T = gas temperature
C pump is computed by a lookup table and multiplied by the speed and pressure to form the initial flow
estimate. During transients, the throttle rate, with the derivative approximated by a high-pass filter,
corrects the air flow for filling dynamics. The control algorithm provides additional correction according
to Eqs. 6.2:
0.5, EGO ≤ 0.5

e0 = 
−0.5, EGO > 0.5
e1 = K i ( N , P m )e 0 Equation 6.2
e , LOW mode with valid EGO signal
e˙ 2 =  1
0 , RICH, DISABLE or EGO warmup
e2 = closed–loop correction
e 0, e1, e2 = intermediate error signals
The nonlinear oxygen sensor, modeled with a hyperbolic tangent in the engine gas Mixing and
Combustion subsystem, provides a meaningful signal when in the vicinity of 0.5 volt. The raw error in
the feedback loop is thus detected with a switching threshold, as indicated in Equation 6.2. If the value is
low (the mixture is lean), the original air estimate is too small and needs to be increased. Conversely,
when the oxygen sensor output is high, the air estimate is too large and needs to be decreased. Integral
control is utilized so that the correction term achieves a level that brings about zero steady-state error in
the mixture ratio.
The normal closed-loop operation mode, LOW, adjusts the integrator dynamically to minimize the error.
The integration is performed in discrete time, with updates every 10 milliseconds. When operating open-
loop however, in the RICH or O2 failure modes, the feedback error is ignored and the integrator is held.
This gives the best correction based on the most recent valid feedback.
Fuel Calculation
The Fuel Calculation subsystem (Figure 6.7) sets the injector signal to match the given airflow
calculation and fault status. The first input is the computed airflow estimation. This is multiplied with
the target fuel/air ratio to get the commanded fuel rate. Normally the target is stoichiometric, 1/14.6.

When a sensor fault occurs, the Stateflow control logic sets the mode input to a value of 2 or 3 (RICH or
DISABLED) so that the mixture is either slightly rich of stoichiometric or is shut down completely.
Fuel Rate Calculation
1
est.
air
flow feedforward fuel rate
1/14.6
F/A Norm
1/(14.6*0.8)
F/A Rich mode
0
mode Shutdown fuel rate 1
4 fuel
limit rate
3 Failures output
Failures
2 feedback correction
correction
Switchable
Compensation
Figure 6.7: Fuel Calculation Subsystem
The Fuel Calculation subsystem (Figure 6.7) employs adjustable compensation (Figure 6.8) in order to
achieve different purposes in different modes. In normal operation, phase lead compensation of the
feedback correction signal adds to the closed-loop stability margin. In RICH mode and during EGO
sensor failure (open loop), however, the composite fuel signal is low-pass filtered to attenuate noise
introduced in the estimation process. The end result is a signal representing the fuel flow rate which, in
an actual system, would be translated to injector pulse times.

2
max
mode
==
2
LOW
3
Failures 8.7696z - 8.5104
O2 fail
z - 0.74082
(warmup)
LOW Mode
1
feedforward
fuel rate
== 1
Fuel Rate
RICH
4
0.25918
feedback
correction z - 0.74082
RICH Mode
Shutoff
Mode
0
Loop Compensation and Filtering

Figure 6.8: Switchable c ompensation
Results and Simulation results are shown in Figure 6.9 and Figure 6.10. The simulation is run with a throttle input
Conclusions that ramps from 10 to 20 degrees over a period of two seconds, then back to 10 degrees over the next two
seconds. This cycle repeats continuously while the engine is held at a constant speed so that the user can
experiment with different fault conditions and failure modes. To simulate a sensor failure, double-click
on its associated switch (see Figure 6.1). Repeat this operation to toggle the switch back for normal
operation.
Figure 6.9 compares the fuel flow rate under fault-free conditions (baseline) with the rate applied in the
presence of a single failure in each sensor individually. Note, in each case, the nonlinear relationship
between fuel flow and the triangular throttle command (shown qualitatively on its Simulink icon). In
the baseline case, the fuel rate is regulated tightly, exhibiting a small ripple due to the switching nature of
the EGO sensor’s input circuitry. In the other four cases the system operates open loop. The control
strategy is proven effective in maintaining the correct fuel profile in the single-failure mode. In each of
the fault conditions, the fuel rate is essentially 125% of the baseline flow, fulfilling the design objective of
80% rich.
Figure 6.10 plots the corresponding air/fuel ratio for each case. The baseline plot shows the effects of
closed-loop operation. The mixture ratio is regulated very tightly to the stoichiometric objective of 14.6.
The rich mixture ratio is shown in the bottom four plots of Figure 6.10. Although they are not tightly
regulated, as in the closed-loop case, they approximate the objective of air/fuel = 0.8(14.6) = 11.7.

FaultTolerant Fuel Control System: Fuel Rate
g/sec
1
baseline
0
g/sec
1
throttle sensor failed
0
2
g/sec
1
speed sensor failed
0
2
g/sec
1
EGO sensor failed
0
2
g/sec
1
MAP sensor failed
0
0 1 2 3 4 5 6 7 8
time (sec)
Figure 6.9: Comparative results for simulated fuel rate
FaultTolerant Fuel Control System: Air/Fuel Ratio
15
baseline
10
15
throttle sensor failed
10
15
speed sensor failed
10
15
EGO sensor failed
10
15
MAP sensor failed
10
0 1 2 3 4 5 6 7 8
time (sec)
Figure 6.10: Comparative results for simulated a ir/fuel ratio

The transient behavior of the system is shown in Figure 6.11. With a constant 12° throttle angle and the
system in steady-state, a throttle failure is introduced at t = 2 and corrected at t = 5. At the onset of the
failure, the fuel rate increases immediately. The effects are seen at the exhaust as the rich ratio propagates
through the system. The steady-state condition is then quickly recovered when closed-loop operation is
restored.
Response to Throttle Failure

15
14
air/fuel ratio
13
12
11
10
2.5
2
fuel rate (g/sec)
1.5
0.5
0
0 1 2 3 4 5 6 7 8
time (sec)
Figure 6.11: Transient response to fault d etection
During simulation, this behavior can also be observed from the Stateflow perspective. By enabling
animation in the Stateflow debugger, the state transitions are highlighted in Figure 6.3 as the various
states are activated. The sequence of activation is indicated by changing colors. This closely coupled
synergy between Stateflow and Simulink fosters the modeling and development of complete control
systems. An engineer’s concepts can develop in a natural and structured fashion with immediate visual
feedback reinforcing each step.

VII. A UTOMATIC T RANSMISSION C ONTROL
Summary In this example, Simulink is used to model an automotive drivetrain. Stateflow enhances the Simulink
model with its representation of the transmission control logic. Simulink provides a powerful
environment for the modeling and simulation of dynamic systems and processes. In many systems,
though, supervisory functions like changing modes or invoking new gain schedules must respond to
events that may occur and conditions that develop over time. As a result, the environment requires a
language capable of managing these multiple modes and developing conditions. In the following
example, Stateflow demonstrates its strength in this capacity by performing the function of gear selection
in an automatic transmission. This function is combined with the drivetrain dynamics in a natural and
intuitive manner by incorporating a Stateflow block in the Simulink block diagram.
Engine Torque Gearset and Vehicle mph

Throttle
Converter Shift Mechanism Dynamics
Transmission
Transmission brake
Control Unit
Figure 7.1: Drivetrain system block diagram
Figure 7.1 shows the power flow in a typical automotive drivetrain. Nonlinear ordinary differential
equations model the engine, four-speed automatic transmission, and vehicle. The model directly
implements these as modular Simulink subsystems. On the other hand, the logic and decisions made in
the transmission controller (TCU) do not lend themselves to well-formulated differential or difference
equations; these are better suited to a Stateflow representation. Stateflow monitors the events which
correspond to important relationships within the system and takes the appropriate action as they occur.
Analysis and The engine receives input in the form of the throttle opening, as commanded by the driver. It is connected
Physics to the impeller of the torque converter which couples it to the transmission
I ei N˙ e = T e − Ti
N e = engine speed
. I ei = engine + impeller moment of inertia Equation 7.1
T e = f1(throttle, N e ) = engine torque
Ti = impeller torque
The input-output characteristics of the torque converter can be expressed as functions of the engine speed
and the turbine speed. In this example, the direction of power flow is always assumed to be from impeller
to turbine.

Ti = (N e / K )2
K = f 2 (N in / N e )
= capacity or K–factor
Equation 7.2
N in = turbine (torque converter output) speed
= transmission input speed
Tt = RTQTi
= turbine torque
RTQ = torque ratio
= f 3 ( N in / N e )
The transmission model is expressed as static gear ratios, assuming small shift times.
RTR = f 4 ( gear )
T out = RTRTin
N in = RTR N out
Equation 7.3
Tin ,T out = transmission input and output torque
N in , N out = transmission input and output speed
RTR = transmission ratio
The final drive, inertia, and a dynamically varying load constitute the vehicle dynamics.
I v N˙ w = R fd (T out − T load )
I v = vehicle inertia
N w = wheel speed
Equation 7.4
R fd = final drive ratio
T load = load torque
= f 5(Nw )
The load torque includes both the road load and brake torque. The road load is the sum of frictional and
aerodynamic losses.
T load = sgn( mph )( R load 0 + R load 2 mph 2 + T brake )

T load = load torque
R load 0 , R load 2 = friction and aerodynamic drag coefficients Equation. 7.5
T brake = brake torque
mph = vehicle linear velocity
The model programs the shift points for the transmission according to a schedule, such as is shown in
Figure 7.2. For a given throttle in a given gear, there is a unique vehicle speed at which an upshift takes
place; the simulation operates similarly for a downshift.

transmission shift points
100
90 upshift
downshift
80
70 4-3
vehicle speed (mph)
60
3-4
50
40
3-2
2-3
30
20
1-2 2-1
10
0
0 10 20 30 40 50 60 70 80 90 100
throttle (%)
Figure 7.2: Shift Schedule
Modeling The Simulink model (sf_car.mdl) is composed of modules which represent the engine, transmission
and vehicle, with an additional shift logic block to control the transmission ratio. User inputs to the
model are in the form of throttle (%) and brake torque (ft-lb). The diagram in Figure 7.3 shows the
overall model.
sf_car.mdl impeller torque

[0 0;100 0;101 400;200 400] brake Choose Start from
schedule
Ti the Simulation menu
to run the simulation.
Ne Ne
engine RPM Ti
9 40;15 100;10 throttle
throttle Engine gear

schedule
Tout
Nout output torque
throttle Vehicle vehicle
gear transmission speed
vehicle_speed transmission speed vehicle mph

(yellow)
shift_logic & throttle %
Mu
Figure 7.3: Overall simulation model
The Engine subsystem consists of a two-dimensional table that interpolates engine torque vs. throttle and
engine speed. In accordance with Equation 7.1, the model subtracts the impeller torque, divides the
difference by the inertia and then numerically integrates the quotient to compute the engine speed.
Figure 7.4 shows the composite engine subsystem.

1
Ti
2 1
1/Iei 1
throttle s
Ne
speed
engine engine + impeller
torque inertia
Figure 7.4: Engine subsystem
The torque converter and the block which implements the various gear ratios make up the transmission
subsystem, as shown in Figure 7.5.
engine speed impeller torque

1 1
Ne Ti
turbine torque
Tin
Tout 2
output torque
Torque gear Tout
Converter Nin
Nout
2
gear transmission
ratio
transmission output speed
3
Nout
turbine speed
Figure 7.5: Transmission subsystem
The torque converter is a masked subsystem, under which the model computes the relationships of
Equation 7.2. The parameters entered into the subsystem are a vector of speed ratios (Nin/Ne ) and vectors
of K-factor (f2) and torque ratio (f3) corresponding to the speed ratio data. Figure 7.6 shows the subsystem
implementation.
1
Ti
Nin
u^2
2
2
impeller Tt
1 turbine
Ne K quotient
speed
ratio factor
Torque
ratio
TORQUE CONVERTER
Figure 7.6: Torque c onverter subsystem
The transmission ratio block determines the ratio RTR(gear), shown in Table 7.1 and computes the
transmission output torque and input speed, as indicated in Equation 7.3. The ratios used progress from
low to another underdrive ratio, one-to-one and overdrive.

gear RTR
1 2.393
2 1.450
3 1.000
4 0.677
Table 7.1: Transmission Gear Ratios
Figure 7.7 shows the block diagram for the subsystem that realizes this ratio in torque and speed.
1
1
Tin
Tout
Product
2
gear
LookUp
Table
2 3
Nin Nout
Product1
Figure 7.7: Transmission gear ratio subsystem
The Stateflow block labeled shift_logic implements gear selection for the transmission. The Stateflow
Explorer is utilized to define the inputs as throttle and vehicle speed and the output as the desired gear
number. Two dashed AND states keep track of the gear state and the state of the gear selection process. The
overall chart is executed as a discrete-time system, sampled every 40 mSec. The Stateflow diagram shown
in Figure 7.8 illustrates the functionality of the block.
The shift logic behavior, explained in the following, can be observed during simulation by enabling
animation in the Stateflow debugger. The selection_state, which is always active, begins by
performing the computations indicated in its during function. The model computes the upshift and
downshift speed thresholds as a function of the instantaneous values of gear and throttle (see Figure 7.2).
While in steady_state, the model compares these values to the present vehicle speed to determine if a
shift is required. If so, it enters one of the confirm states (upshift_confirm or downshift_confirm),
which records the time of entry.
If the vehicle speed no longer satisfies the shift condition, while in the confirm state, the model ignores the
shift and it transitions back to steady_state. This prevents extraneous shifts due to noise conditions. If
the shift condition remains valid for a duration of Tconfirm, the model transitions through the lower
junction and, depending on the current gear, it broadcasts one of the shift events. Subsequently, the model
again activates steady_state after a transition through one of the central junctions. The shift event,
which is broadcast to the gear_selection state, activates a transition to the appropriate new gear.

shift_logic
gear_selection/
gear_state
UPSHIFT34
UPSHIFT12 UPSHIFT23
first/ second/ third/ fourth/
entry: gear = 1; entry: gear = 2; entry: gear = 3; entry: gear = 4;
DOWNSHIFT21
DOWNSHIFT32 DOWNSHIFT43
selection_state/
during: down_threshold = ml(’interp2([1:4],downth, downtab, %g, %g)’, gear, throttle);...
up_threshold = ml(’interp2([1:4],upth, uptab, %g, %g)’, gear, throttle);
[vehicle_speed > up_threshold]

downshift_confirm/ steady_state upshift_confirm/
entry: tdn = t; entry: tup = t;
[vehicle_speed < down_threshold]
[vehicle_speed > down_threshold]

[vehicle_speed < up_threshold]
[t-tdn >= Tconfirm & ... [t-tup >= Tconfirm & ...
vehicle_speed <= down_threshold] vehicle_speed >= up_threshold]
[gear == 3]/UPSHIFT34
[gear ==4]/DOWNSHIFT43
[gear == 2]/UPSHIFT23
[gear == 3]/DOWNSHIFT32
[gear == 2]/DOWNSHIFT21 [gear == 1]/UPSHIFT12
Printed 11Sep1997 16:52:48
Figure 7.8: Stateflow diagram of the transmission shift logic
For example, if the vehicle is moving along in second gear with 25% throttle, the state second is active
within gear_state, and steady_state is active in the selection_state. The during function of the latter
finds that an upshift should take place when the vehicle exceeds 30 mph. At the moment this becomes
true, the model enters the upshift_confirm state and sets the local variable tup to the current time by its
entry action. While in this state, if the vehicle speed remains above 30 mph until the elapsed time (t -
tup) reaches Tconfirm (0.1 Sec), the model satisfies the transition condition leading down to the lower
right junction. This also satisfies the condition [gear == 2] on the transition leading from here to
steady_state, so the model now takes the overall transition from upshift_confirm to steady_state
and broadcasts the event UPSHIFT23 as a transition action. Consequently, the transition from second to
third is taken in gear_state which completes the shift logic.
The vehicle dynamics (Figure 7.9) use the net torque to compute the acceleration and integrate it to
compute the vehicle speed, per Equation 7.4 and Equation 7.5. In this example, we again use a masked
subsystem for the vehicle submodel. The parameters entered in the mask menu are the final drive ratio,
the polynomial coefficients for drag friction and aerodynamic drag, the wheel radius, vehicle inertia, and
initial transmission output speed.

2 Rfd
1/Iv 1/s Rfd 2
transmission
output torque final wheel transmission
drive vehicle final drive
inertia speed ratio output speed
ratio
1
vehicle
linear speed mph speed
2*pi*Rw 60/5280 f(u)

road load
1
brake torque signed
load
Figure 7.9: Vehicle dynamics subsystem
Results The engine torque map, torque converter characteristics, and road load data used in the simulations are
shown in the three plots which follow (Figure 7.10, Figure 7.11, and Figure 7.12).
Engine Torque Map f1

350
100%
300 70%
60%
50%
250
40%
200
engine torque (ft-lb)
30%
150
100 20% throttle
50
0% throttle
-50
-100
500 1000 1500 2000 2500 3000 3500 4000 4500 5000
engine speed (RPM)
Figure 7.10: Engine map

Torque Converter Characteristics, f , f
2 3
260
240
220
2.0
200
K-factor (RPM/sqrt(ft-lb))
Torque Ratio
180
160
140 K-factor
120
1.0
100
80
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
speed ratio, N /N
in e
Figure 7.11: Torque c onverter c haracteristics
Vehicle Road Load, f5

250
200
torque (ftlb)
150
100
50
80
60
power (HP)
40
20
0
0 10 20 30 40 50 60 70 80 90 100
vehicle speed (mph)
Figure 7.12: Vehicle road load or d rag torque

The first simulation uses the following throttle schedule:
 0 60 
 
14.9 40 
throttle =  15 100 
 
 100 0 
 200 0 
 
The first column corresponds to time; the second column corresponds to throttle opening in percent. In
this case we do not apply the brake (0 ft-lb). The vehicle speed starts at zero and the engine at 1000 RPM.
Figure 7.13 shows the plot for the baseline results, using the default parameters. As the driver steps to
60% throttle at t = 0, the engine immediately responds by more than doubling its speed. This brings
about a low speed ratio across the torque converter and, hence, a large torque ratio (see Figs. 7.6 and 7.11).
The vehicle accelerates quickly (no tire slip is modeled) and both the engine and the vehicle gain speed
until about t = 2, at which time a 1-2 upshift occurs. The engine speed characteristically drops abruptly,
then resumes its acceleration. The 2-3 and 3-4 upshifts take place at about four and eight seconds,
respectively. Notice that the vehicle speed remains much smoother due to its large inertia.
Baseline Automatic Transmission Simulation

5000
engine speed (RPM)
4000
3000
2000
1000
140
vehicle speed (mph) & throttle (%)
120
100
80 vehicle mph
60
throttle %
40
20
0
0 5 10 15 20 25 30
time (sec)
Figure 7.13: Initial simulation time history
At t = 15, the driver steps the throttle to 100% as might be typical of a passing maneuver. The
transmission downshifts to third gear and the engine jumps from about 2600 to about 3700 RPM. The
engine torque thus increases somewhat, as well as the mechanical advantage of the transmission. With

continued heavy throttle, the vehicle accelerates to about 100 mph and then shifts into overdrive at about t
= 21. The vehicle cruises along in fourth gear for the remainder of the simulation.
Figure 7.14 shows the results of a second simulation. The behavior for the first fifteen seconds is the same
as above, but the throttle subsequently drops to about 5% at 40 seconds. This is followed by a step in brake
torque at t = 50. Again, the large vehicle inertia dominates the dynamics as it eventually slows down to a
crawl. The engine speed downshifts occur at about 72, 80 and 90 seconds, ending in first gear.
Acceleration and Braking Maneuvers

6000
engine speed (RPM)
4000
2000
150
vehicle speed (mph)
100
50
0
throttle (%) & torque (ftlb)
300
brake
200
throttle
100
0
0 10 20 30 40 50 60 70 80 90 100
time (sec)
Figure 7.14: Vehicle simulation with a cceleration and b raking
Conclusions We can easily enhance this basic system in a modular manner, for example, by replacing the engine or
transmission with a more complex model. We can thus build up large systems within this structure via
step-wise refinement. The seamless integration of Stateflow control logic with Simulink signal
processing enables the construction of a model which is both efficient and visually intuitive.

VIII. E LECTROHYDRAULIC S ERVO CONTROL
Summary In this example we develop a Simulink model for a hydraulic servomechanism controlled by a pulse-
width modulated (PWM) solenoid. This might represent a motion control system in an industrial or
manufacturing setting, or a subsystem that controls the position of a valve in an automotive or aerospace
application. Nonlinear differential equations are used to model the magnetic, hydraulic and mechanical
components; discrete-time difference equations represent the controller. A behavioral model in Stateflow
implements the electronic circuit which generates the PWM waveforms and regulates the solenoid
current. Although a detailed power electronic model could be developed in Simulink, the Stateflow
description provides the required functionality and speeds development.
Figure 8.1 shows a hydraulic schematic for the mechanism. The objective of the system is to position the
load xp so that it follows commands issued in the form of a time-varying set point rset . An electronic
controller compares these commands to feedback measurements of xp and generates a PWM control
signal at a rate of 50 Hz. The PWM duty cycle is the percentage of the 20 millisecond period for which the
valve directly supplies oil to the control pressure developed in the cylinder behind the piston, pc. For the
remainder of the period, the valve vents pc to exhaust. The composite flow qnet thus controls pc which
develops an actuating force against the piston. This forces the spring-loaded piston to its position xp such
that it follows the reference trajectory rset .
Cylinder and Piston

qnet
pc
Ap
Ksp
PWM
Ps Solenoid xp
Valve
Figure 8.1: Solenoid v alve and h ydraulic actuator
We chose PWM control to regulate the net valve flow with its on/off duty ratio rather than relying on the
strict mechanical tolerances of a continuous valve. This sequentially turns the solenoid completely on or
off , rather than attempting to control it to a precise intermediate position. The tradeoff is that a
disturbance is introduced into the system at the PWM frequency. As a result, we must take care that this
is adequately attenuated by the low-pass response of the mechanical system. We can evaluate this
requirement by constructing a simulation at the design phase rather than waiting for experimental parts.
Analysi s and PWM S OLENOID
Physics The model of the solenoid-controlled PWM valve includes three parts:
Magnetic circuit
Armature motion
Valve flows

Figure 8.2 shows a cross-sectional view of a typical solenoid valve of this type. The enclosure, armature
and pole piece are steel, and the coil is wound around the armature/pole axis. With no current, the
internal spring forces the armature and ball to the right against the hydraulic force. This blocks the
supply pressure Ps and opens a path from control pressure to exhaust. When the solenoid is energized, the
armature and pole come together and the pressure force shuttles the ball to open the supply port and block
the exhaust port.
flux path
supply
pre ssure
pole armature Ps
air gap
exh
control
coil pc pressure
Figure 8.2: Pulse-Width Modulated Solenoid Valve
Consider first the magnetic circuit. Faraday’s law determines the flux. We assume that fringing and
leakage flux are negligible, as are eddy currents.
φ˙ = (vsol − iR )/ N ,
φ = flux
vsol = solenoid voltage
i = current
R = winding resistance
N = number of turns
The magnetomotive force required to develop this flux is broken up into components for the steel and the
air gap. Although the majority of the circuit’s reluctance is concentrated at the air gap, the nonlinear
properties of the steel components, such as saturation and hysteresis, can limit performance.

MMF = MMFair + MMFsteel
MMFair = H air g
MMFsteel = Hsteel Lsteel
MMF = magnetomotive force Equation 8.1
H = magnetic field intensity
g = length of air gap
Lsteel = magnetic circuit length in steel
Within the steel, the flux density, B, is a nonlinear function of H, dependent upon the material properties.
We also assume that the area, A, which relates φ and B at the air gap, applies uniformly for the steel path.
B = φ /A
= flux density
= f(Hsteel )
= µ0 H air
A = cross- sectional area at air gap
µ0 = permeability of air
The solenoid force, Fsol , and current that result are:
F sol = 0.5B 2 A / µ 0
i = MMF / N
The armature responds to the solenoid force, as well as the hydraulic and spring forces.
mx˙˙ = Fsol + AoPs − K s x − C v x˙

x = armature position
= g max − g
m = mass
Equation 8.2
Ao = supply orifice area
Ps = supply pressure
K s = return spring rate
C v = damping rate
The net oil flow directed from the valve to the actuator, qnet , is the supply flow less the exhaust flow.
q net = q s − q ex
K A sgn(P − p ) P − p , x > 0

qs =  o o s c s c
0, x =0
K A p , x < balltravel
q ex =  o o c Equation 8.3
0, x = balltravel
p c = control pressure
K o = flow coefficient

Actuator Dynamics
The actuator assembly moves the piston against a spring as a function of the control pressure developed
behind it. Assuming negligible leakage,
β
(
p˙ c = q net − x˙ p Ap ,
V
)
β = fluid bulk modulus
V = x p Ap Equation 8.4
= fluid volume
x p = piston position
Ap = actuator (piston) area
The actuator’s equation of motion, dominated by the relatively large hydraulic and spring forces, is
simply:
M p x˙˙p = p c Ap − K sp x p ,
M p = net actuator mass Equation 8.5
K sp = spring rate
Electronic Controls
We employ a discrete-time PI (proportional + integral) control law to
1. Achieve zero steady-state error to step changes in the position set point, and
2. Compensate for the low-frequency actuator dynamics to improve response speed.
 K 
dutycycle =  K p + I  (rset − x p ) Equation 8.6
 z − 1
The integral term is essential because the null duty cycle, or equilibrium control input, is subject to
uncertainty and will change with the system’s operating point. The proportional part contributes phase
lead at low frequency which is essential for stability.
Equation 8.6 computes the PWM duty cycle as a function of position error. The duty cycle is applied to a
50 Hz pulse train and the power electronics convert the pulse signal to solenoid current. Digital and
analog integrated circuits are available to perform these functions, so we use a behavioral model, rather
than a highly detailed physical model. The behavior is best described in terms of the circuit’s reaction to
the commands it receives and the response of its load. Figure 8.3 shows an idealized example.
At the beginning of each 20 millisecond period, the PWM pulse turns on and must pull the solenoid
armature up against the pole piece to open the valve to supply pressure. Hence, the driver circuit applies
the full supply voltage to achieve the fastest initial rise in current. The solenoid maintains this condition
until the current has risen to the level at which the magnetic and hydraulic forces overcome the spring
and move the armature.

Once the armature has been pulled in, the air gap is very small and somewhat less current is needed to
hold the armature in place. The driver thus regulates the current at a lower level for the remainder of the
“on” portion of the cycle. Typically, a switch-mode regulator controls the “hold” current. This technique
alternatively applies the supply voltage to the solenoid and then allows the field to collapse slowly
(shunted by a freewheel diode). This is significantly more efficient, in terms of power, than linear
regulation.
At the end of each pulse, the armature releases so that the ball returns to its original position and the
valve opens to exhaust. We achieve this by opening the solenoid circuit so that the magnetic field
collapses quickly. Typically, we employ a zener diode to limit the large negative EMF while still allowing
a fast decay. The current then remains at zero for the duration of the “off” time until the next cycle
begins.
In this way, we divide the “on” portion of each pulse into two phases: “pull-in” and “hold.” The “off”
portion is characterized by the initial rapid decay, followed by zero voltage and current. The diagram in
Figure 8.3 illustrates this scenario.
vol tage
Vs Ipull
current
Ihold
20 mSec
pull in hold off
Figure 8.3: Current c ontrol within pulses
Modeling Figure 8.4 shows the top-level system block diagram (sf_electrohydraulic.mdl). The set point block
consists of a signal generator and a step function which are added to give wide flexibility to the user in
specifying the set point. The controller is a straightforward discrete-time subsystem. We implemented the
PWM driver circuit in Stateflow, but it functions just like the other subsystems at the block diagram level.
We refined the solenoid valve into three parts, as described above. The actuator model, consisting of the
cylinder pressurization and piston motion subsystems, completes the overall system model.

Supply cylinder
Pressure control pressure pressurization
Ao piston Mux
Ps Fball
net flow qnet position
xp
stroke armature position set point
xp pc pc
set point set point control (yellow)
rset x Fsol Fsol valve flows pressure xpdot
xpdot and
duty cycle duty_cycle piston position
armature piston motion
position v vsol isol
feedback solenoid motion
command i voltage
Magnetic
Controller PWM_driver_ckt Circuit
solenoid
current (A)
electrohydraulic servomechanism
Figure 8.4: Servo model using Simulink and Stateflow
Controller
The controller samples the position error and generates a new solenoid duty cycle every 20 milliseconds.
The duty cycle consists of a component that is proportional to the error plus a component that is
proportional to the integral of the error. The model realizes the integration in the z domain with feedback
around a 1/z block that places a pole at z = 1. The integral gain, as labeled in the diagram, (Figure 8.5) is
fixed with respect to the proportional gain. An overall loop gain, Ka, adjusts both while keeping their
ratio, hence the transfer function zero, constant. While in the linear operating range:
duty cycle  KI 
= K a 1 + 
error  z − 1 Equation 8.7
 z − (1 − K I ) 
= Ka  
 z −1 
The model limits the computed duty cycle so that it never falls below the minimum time to open the
valve, nor exceeds the time at which the valve remains continuously open. Whenever it reaches either of
these limits, the integrator holds constant (zero input) until the error is of the appropriate sign to pull it
away from the limit.
1
Ka
set point
proportional + integral control
sampled
loop 1
error
gain
0.1332 duty cycle
1 Saturation
2 feedback z
0 integral
gain
dc
nothold
e1
prevent
windup
controller
Figure 8.5: Discrete-time c ontroller subsystem

PWM Driver Circuit
The solenoid driver circuit uses the computed duty cycle to generate the PWM waveform. The solenoid
voltage is applied in order to achieve the desired current, force, and hence, valve flow. We modeled this
with a Stateflow block, which uses duty cycle and solenoid current as inputs and computes voltage as an
output. Figure 8.6 shows the Stateflow diagram for the model.
PWM_driver_ckt
/ton = 0; energize_solenoid
/toff = ton + duty_cycle*Tpwm/100;
pull_in_current/
entry : v = Vs;
[i >= Ipull]
[t >= ton]
solenoid_off/
entry : v = -(i > 0)*Vz; regulate_hold_current
during: v = -(i > 0)*Vz; [t > toff]/ton += Tpwm;
freewheel/
entry: v = -(i > 0)*Vd;
during: v = -(i > 0)*Vd;
[ i <= Ihold-deltai] [i >= Ihold + deltai]
hold/
entry : v = Vs;
Printed 12Sep1997 12:59:58
Figure 8.6: Stateflow diagram for the PWM d river c ircuit
Each PWM cycle begins with the local variable ton equal to the current simulation time. The
unconditional transition which begins the cycle computes toff, the time at which the “on” portion of the
pulse ends.
toff = ton + duty_cycle*Tpwm/100; Equation 8.8
Tpwm is a MATLAB workspace variable representing the pulse period. The system enters the
energize_solenoid state and, by default, the pull_in_current state. As described above, the driver
circuit connects the supply voltage to the output for this phase of the pulse. Once the current reaches
Ipull, the worst-case current required to pull in the armature, it enters the regulate_hold_current
state. A diode in the freewheel state shunts the coil which clamps the solenoid voltage at -Vd. When the
current falls to the hold level, the system alternates between the hold and freewheel states to regulate it to
Ihold ± deltai.
When the time reaches t = toff, it exits the energize_solenoid state, regardless of which of the pull-in,
freewheel or hold states is currently active. This is achieved by drawing the transition directly from the
superstate boundary to the solenoid_off state. The value of ton, the beginning of the next cycle, is updated
at this time. While in the solenoid_off state, the coil connects to the zener voltage, -Vz, until the field
collapses and the current falls to zero, as described in the entry and during actions.

Magnetic Circuit
The model uses the applied voltage and armature position to determine the solenoid force and current.
This requires evaluating Equations 8.1 with Simulink blocks placed in the appropriate configuration. The
state variable is flux, computed by integrating the solenoid EMF. The flux density is calculated by
dividing the flux by the cross-sectional area of the magnetic path. The force Fsol is computed as a gain
times the square of the flux density.
net air gap air mmf overall

gmax
mmf
x gap at x = 0
1
armature position gser
equiv. series gap
2
1
vsol 1/N 1/A 1/mu0
s
R
flux 1/N 2
flux density
isol
current
Lsteel
H vs. B steel mmf
2 0.5*A/mu0
u 1
Fsol
force
magnetic circuit
Figure 8.7: Magnetic subsystem model
The model computes the solenoid current by determining the magnetomotive force. In the air gap, Hair =
B/µ0 is multiplied by the gap length to give MMFair . The gap length is computed by subtracting the
armature position from the maximum gap. A small additional gap is added to model additional air in
the circuit, at the armature bearing surface, for example. The MMF required to produce the flux density in
the steel is computed by putting the material characteristics, H vs. B, in a 2-D lookup table. Since this
curve has significant hysteresis, two curves are placed in the table, one for increasing and one for
decreasing flux. The appropriate curve for H steel = f ( B ) is selected according to the sign of φ̇ . Hsteel Lsteel
is added to MMFair and the sum is divided by N to determine the solenoid current.
Armature Motion
The model solves the equation of motion for the armature directly, as shown in Figure 8.8. The sum and
gains use standard blocks , and the subsystem Double Integrator computes the velocity and position of
the armature based on its acceleration. The position x = 0 corresponds to the maximum air gap, gmax.

1
Fball
2
x 1
Fsol
1/m xdotdot stroke
xdot
mass
double integrator
net force
Cv
Fs0 damping
spring preload
Ks
spring rate
armature motion
Figure 8.8: Mechanical subsystem for the solenoid Armature
In the double integrator subsystem (Figure 8.9), the model limits the position integrator by physical stops
at x = 0 and at x = gmax - gmin (a shim typically limits the minimum gap). When these limits are
reached, it is essential that the velocity becomes zero and remains zero while at the stops. The model
achieves this by feeding the position saturation port back to the velocity reset trigger. In addition, the
derivative input of the velocity integrator switches to zero as long as xdotdot (force/mass) holds the
armature against the stop. The velocity thus remains zero until the force reverses direction.
2
xdot
0
1 1 1
1
s s x
xdotdot Switch
>
0 AND
OR
AND
<
0 <=
0
>=
gmax-gmin
Figure 8.9: Cascaded integrators with c oordinated limit logic
Valve Flows
Simulink models the turbulent flow through the valve orifices with the following subsystem, shown in
Figure 8.10. The inputs pup and pdown are the upstream and downstream pressures and q is the flow
from pup to pdown. The square root of the absolute value of the pressure drop, multiplied by the sign of the
pressure drop and KoAo (defined in Equation 8.3) yields the flow.
pup
1
2 Ko*Ao 1
pdown |u| sqrt
q
Turbulent Orifice Flow
Figure 8.10: Individual o rifice subsystem

The valve model shown in Figure 8.11 uses this subsystem twice, to model the flow from the supply to the
control pressure and to model the flow from the control pressure to the exhaust. The net flow to the
control pressure is the supply flow, when x > 0, and negative one times the exhaust flow, when x <
balltravel. The maximum ball travel is less than the armature travel in order to assure that the ball seats
against the exhaust port when the armature is pulled in. During the brief time in which the ball is at
intermediate positions, with neither port blocked, flow occurs at both orifices.
source flow
Ps pup
q
pdown
armature position
2
>
0 source
open
0
1
net
>= flow
balltravel exhaust
closed
1 pup
control pressure q
0 pdown
exhaust flow
poppet valve flows
Figure 8.11: Overall v alve flow subsystem
Cylinder Pressurization
The model for cylinder pressurization is a direct realization of Equation 8.4 in the actuator dynamics
section (see Figure 8.12). Oil volume is the product of the piston position and its cross-sectional area. The
division operator uses a function block within a masked subsystem and the other blocks are standard
gains, a sum, and an integrator.
bulk
modulus
1
beta
qnet 1
1
effective s
flow pc
pressure
3 Ap
xpdot
dVolume/dt
2 Ap
xp
Volume
cylinder pressurization
Figure 8.12: Hydraulic c ylinder subsystem

Piston Motion
Equation 8.5 is the differential equation for piston motion in the previously stated actuator dynamics
section. The Simulink implementation is straightforward, as shown in Figure 8.13. We again use the
double integrator subsystem, described above in the armature motion section, to insure that zero velocity
is indicated when the actuator is being held against its physical stops.
spring
Ksp
xp
x 1
Ap 1/Mp xdotdot
1
pc
xdot 2
area mass xpdot
double integrator
piston motion
Figure 8.13: Actuator mechanical subsystem
Results Figure 8.14 below shows the set point and piston position for a baseline simulation. During the first 0.1
second, and again from 1.0 to 1.1 seconds, the output is slew rate limited by the maximum flow available
to the actuator. At other times, the 3 Hz sinusoid is tracked closely. Although the solenoid goes through a
complete on/off cycle each PWM period, the 50 Hz dither superimposed on the actuator position is
relatively small.

Electrohydraulic Servo Control
0.018
0.016
set point
0.014
0.012
actuator position (m) position

0.01
0.008
0.006
0.004
0.002
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
time (sec)
Figure 8.14: Simulated p iston motion
Figure 8.15 below depicts the solenoid current control, under the authority of the Stateflow model. The
diagram shows two cycles with about 47% and about 55% duty cycle, respectively. During the pull-in
phase, as the flux builds and the current approaches its 2.5A target, the current drops abruptly at about t =
2 milliseconds. This is the instant at which the armature is pulled in. This pull-in generates so much
back EMF that the current drops appreciably. The notch in current is so distinct that it is often used in the
laboratory to measure solenoid response time.
When the current reaches the conservative 2.5A target, more than enough to achieve armature pull-in,
the solenoid enters the hold phase of its energized state. The model regulates the average current to 1A by
chopping the voltage as described in the Stateflow diagram. The chopping takes place at a rate somewhat
higher than 1 kHz in order to regulate the current within ±0.1 A. The freewheel state uses a value of Vd =
0.5 V to slow the decay of energy in the magnetic field. When the completion of each energized state turns
off the solenoid, the negative voltage is limited by Vz = 50 V. The model achieves a rapid decay in current
without subjecting the semiconductor devices to extreme voltages.

Electrohydraulic Servo Control
2.5
PWM Driver Current Waveform
solenoid current (A)

1.5
0.5
0
0.94 0.945 0.95 0.955 0.96 0.965 0.97 0.975
time (sec)
Figure 8.15: Simulated solenoid c urrent
Conclusions Simulink and Stateflow combine to provide a powerful modeling environment for dynamic systems. In
this case, Simulink enables the direct construction of block diagram subsystems which represent the
nonlinear differential equations of the physical system and the difference equations of its discrete-time
controller. A Stateflow model captures the behavior of the electronic PWM driver circuit, without
resorting to the complexity of a detailed circuit model. The clear and natural logic of Stateflow facilitates
rapid model development and debugging. The overall model develops in a structured, hierarchical
manner, amenable to careful documentation

IX. MODELING S TICK -S LIP F RICTION
Summary The model in this example consists of a block sliding along a surface and compressing a spring under the
influence of a user-designated input force. In the absence of friction, this behaves like a classical spring-
mass system with the steady-state block position proportional to the applied force. When friction is taken
into account, the model becomes considerably more complicated. Friction between the block and surface
tends to resist motion; however, the friction force changes with velocity and tends to be greatest when
stationary. This results in motion which alternately “sticks” and “slips” as the overall force balance
requires. This “stiction” phenomenon is common in many mechanical systems.
In this simulation, Stateflow is used to represent some of the physical states of the system. As noted above,
the friction force between two surfaces is intrinsically tied to their instantaneous relative velocity. The
continuous trajectory of velocity and position is subject to abrupt changes in acceleration, however,
corresponding to transitions between the discrete states of “stuck” and “sliding.” Simulink provides a
powerful tool for modeling the continuous dynamics, and Stateflow is a natural and intuitive setting for
modeling the discrete physical states.
Analysis and The diagram in Figure 9.1 shows the mechanical system.
Physics
Fin K
M
Figure 9.1: Spring-Mass-Friction System
The basic equation of motion for the block is:
Mx˙˙ = Fin − Fspring − F friction

M = block mass
Equation 9.1
x˙˙ = acceleration
Fin = input force
The model for the linear spring (of negligible mass) is:

F spring = Kx
K = spring rate Equation 9.2
x = position
The friction force is more complex, however.
F friction = 
()
sgn x˙ µF n , Fstationary > µF n
Fstationary , otherwise, at x˙ = 0
x˙ = velocity
()
µ = µ x˙ = coefficient of friction Equation 9.3
F n = normal force
Fstationary = instantaneous force such that x˙ = 0
In many applications the friction capacity is described by its static and kinetic magnitudes. This
approach is used in the present model, also assuming a constant normal force.
 µ static Fn = F static , x˙ = 0
µFn =  Equation 9.4
µkinetic Fn = F sliding , x˙ ≠ 0
The following logic determines Fstationary . Whenever the velocity is nonzero, an impulsive force would be
needed to make it zero instantaneously. This always exceeds the capacity, Fsliding , so the latter magnitude is
used. When the velocity is already zero, however, Fstationary is the force which maintains this condition by
making the acceleration zero.
F stationary = F in − F spring
Equation 9.5
= F sum
The friction force can thus be expressed as
sgn(x˙ )F x˙ ≠ 0
sliding ,


F friction = F sum , x˙ = 0 , F sum < F static Equation 9.6

sgn(F sum )F static , x˙ = 0 , F sum ≥ F static
Modeling Figure 9.2 illustrates a Simulink model that demonstrates this behavior (sf_stickslip.mdl). The two
main components are the block labeled Mechanical Motion and the block labeled state_logic . The former is
composed of conventional hierarchical Simulink subsystems and the latter is implemented in Stateflow.
Simulink is ideal for solving ordinary differential equations and the associated linear and nonlinear
signal flow calculations. Stateflow demonstrates its power in its ability to recognize system events which
require changes in the mode of operation.

sf_stickslip.mdl
Mu force and
position
The input force linearly compresses vs. time
the spring, but friction resists this Mux
movement. The magnitude of position
friction depends on the state of motion. vs.
force
position
Fin zero threshold
Input velocity novelocity

Force
stuck
stuck
Fsum Fsum
mechanical state_logic
motion
0 With the default parameters, the natural

edit frequency is much higher than that of the
excitation force. For contrast, change
parameters the parameter values to:
M = 0.1 kg and Fsliding = 0.1 N.
t
Timekeeping
Stickslip Friction Simulation

To run, choose Start from the Simulation menu.
Figure 9.2: Simulink b lock d iagram
As described in Equation 9.2, the model implements the fundamental equation of motion, in the
Mechanical Motion block. Double-clicking on this block shows the underlying subsystem, pictured in
Figure 9.3. The sum of the forces divided by the mass determines the block acceleration. The acceleration
is integrated twice to compute the velocity and position.
Fsum
3
2
1
velocity
Fin 1/M 1 1
1
Fsum s s
position
xdot Ffriction
stuck
friction force
2
stuck
spring
Mechanical Motion
Figure 9.3: Mechanical subsystem
The friction force subsystem shown in Figure 9.4 performs a number of nonlinear operations on the
signals in order to model the relationships of Equation 9.6. Standard Simulink blocks implement the

functions of absolute value, sign, minimum, and product. The switch block selects the appropriate value
for the friction force, under control of the signal labeled “stuck.” This is the output of the Stateflow
control_logic block. Note that it is also used in the Mechanical Motion block as a reset input to the first
integrator. This is to ensure that, at the onset of the zero velocity mode, any infinitesimal value is cleared
from the velocity state.
1
Sign1
Fsum |u|
Abs fstatic
min
Fstatic
MinMax
3 1
stuck Ffriction
Fsliding
fsliding
2
xdot
friction forces
Figure 9.4: Friction Subsystem
The Stateflow diagram in Figure 9.5 below describes the behavior of the system states. The block input
signals are Fsum and novelocity. Fsum as defined in Figure 9.4 and novelocity is a binary signal that
becomes one the instant the velocity crosses through zero. The Stateflow block output is the control signal
stuck, as described above. The parameter Fstatic is taken from the MATLAB workspace.
state_logic
state_of_motion
[fabs(Fsum) > Fstatic]
stuck/ sliding/
entry: stuck = 1; entry: stuck = 0;
[novelocity & (fabs(Fsum)<= Fstatic)]
Printed 14 Nov 1997 14:52:39
Figure 9.5: Stateflow diagram for stick-slip motion

Two mutually exclusive (OR) states are used to represent the conditions of stuck and sliding. The
system, assumed to initially be at rest, first enters the stuck state via a default transition. The transition
out of the stuck state, the upper arc from left to right, is activated when the external forces exceed the static
friction. The sliding state remains active as long as the block is moving, its velocity is nonzero. When
the velocity reaches zero, its direction of travel will reverse, unless the net external force is lower in
magnitude than the static friction force. Hence, a transition from the sliding state to the stuck state is
made when the logical and of these two conditions, zero velocity and force less than static friction, is
satisfied.
The output signal of the machine, stuck, is a binary representation of the state. The entry function of the
stuck state sets it to one, and the entry function of the sliding state clears it to zero. It can then be used as
a control signal in the other Simulink blocks, as indicated above. The ultimate value of the state machine
in this example is to model the state transitions of the system in accordance with physical laws while
using the appropriate frictional force in the acceleration computations at each instant.
Results The default parameters used in this model are:

M = 0.001 kg
K = 1 N/m
F static = 1 N
F sliding = 1 N
The input force ramps linearly from zero to 5 N and back to zero, with a period of 5 seconds. Two
noteworthy characteristics of the parameters are:
1. The natural frequency of the system
ω n = K / M = 31.6 rad/sec
is much higher than that of the excitation (2π/10 rad/sec). As is typical in the laboratory, the model
employs a very low excitation frequency in order to determine the static response characteristics of the
system.
2. The static and kinetic friction magnitudes are equal.
Figure 9.6 and Figure 9.7 show plots of the simulation results. Figure 9.6 shows the time histories of the
input force and the resulting position. The input force must exceed that of the static friction in order to
begin motion at t = 1. From 1 < t < 5, the position tracks the spring force less the kinetic friction force, with
small oscillations showing changes in velocity at the natural frequency ω n . The input force then begins
to decrease. The mass immediately comes to a halt and sticks until t = 7, when the net force again exceeds
the static friction force, now in the reverse direction.

Baseline Simulation
5
Force
4.5
3.5
Force (N) and position (m)

3
2.5
Position
2
1.5
0.5
0
0 2 4 6 8 10 12 14 16 18 20
time (sec)
Figure 9.6: Simulated time histories
Similar behavior follows as the input force decreases to zero and repeats another cycle. Figure 9.7 shows
the input-output characteristics of the system, position vs. force. The plot forms a hysteresis loop as the
position lags the force. The static input-output characteristics cannot be represented by a single-valued
function because the system has memory.
Friction Characteristics
5
3
position (m)
-1
-1 0 1 2 3 4 5 6
force (N)
Figure 9.7: Simulated hysteresis loop

The continuous states, position and velocity, have memory in the sense that the model stores the energy
according to their magnitudes. The potential energy in the spring is proportional to the square of position
and the kinetic energy of the mass is proportional to the square of its velocity. They represent memory in
that they cannot change instantaneously if power flow is assumed to be finite.
The static behavior relies, not only on the spring rate (force vs. position), but also on the position and
velocity at the instant of the last state transition. While in the stuck state, the position will remain
constant at the point where it entered the state. In the sliding state, the position depends on the spring
characteristics, the direction of the velocity, and the position of the mass at the moment it began to slide.
This behavior is captured in a natural and intuitive way in the Stateflow model.
We can illustrate the dynamic behavior of “stiction,” or stick-slip friction, even more dramatically by
changing the system parameters as follows.
M = 0.1 kg
F sliding = 0.1 N
Stiction Nonlinearity
6
4
force (N) and position (m)
Force
2
Position
0
-1
0 2 4 6 8 10 12 14 16 18 20
time (sec)
Figure 9.8: Stick-slip friction b ehavior
With kinetic friction lower in magnitude than static friction, abrupt discontinuities in acceleration occur
at the stuck-to-sliding and sliding-to-stuck state transitions. As the velocity reaches zero, the
acceleration is often nonzero. If the stuck state is entered, however, the acceleration becomes zero

immediately. This highly nonlinear behavior is typical in many systems, making it difficult to precisely
control position.
Conclusions We can greatly simplify the modeling and simulation task by inserting a Stateflow block into the
mechanical system. Conceptually, Stateflow captures the complex dynamic nonlinear behavior in a
graphical diagram with straightforward, easy-to-read functionality. We can insert this diagram directly
into the Simulink diagram, and the tasks of code generation, compiling, and linking are seamlessly
automated to provide a powerful simulation environment.

I NTERNATIONAL D ISTRIBUTORS AND R ESELLERS
Australia India, Sri Lanka New Zealand South Africa
CEANET Pty., Ltd. Cranes Software International (P) Ltd. Hoare Research Software Opti-Num Solutions
Tel: +61 (0)2-9922-6311 Tel: +91 (0)80-5549-338 Tel: +64-7-839-9102 Tel: +27-11-325-6238
Fax: +61 (0)2-9922-5118 Fax: +91 (0)80-5546-299 Fax: +64-7-839-9103 Fax: +27-11-325-6239
E-mail: belinda@ceanet.com.au E-mail: matlab@CRANES.XEEBLR. E-mail: info@hrs.co.nz E-mail: info@optinum.co.za
Web: www.ceanet.com.au xeemail.ems.vsnl.net.in Web: www.hrs.co.nz Web: www.optinum.co.za
Brisbane office: Authorized Distributor Authorized Reseller Authorized Reseller
Tel: +61 (0)7-3369-4499
Israel The Nordic Countries and Spain, Portugal
Fax: +61 (0)7-3369-4469
Omikron Delta (1927) Ltd. Baltic States Addlink Software Cientifico
E-mail: ken@ceanet.com.au
Tel: +972 (0)3-561-5151 Computer Solutions Europe AB Tel: +34 (9)3 415-49-04
Authorized Distributors
Fax: +972 (0)3-561-2962 Tel: +46 (0)8-15-30-22 Fax: +34 (9)3 415-72-68
Brazil E-mail: info@omikron.co.il Fax: +46 (0)8-15-76-35 E-mail: info@addlink.es
OpenCadd Computacao Grafica Web: www.omikron.co.il E-mail: info@comsol.se Web: www.addlink.es
Tel: +55-11-816-3144 Authorized Distributor Web: www.comsol.se Authorized Distributor
Fax: +55-11-816-7864 Søborg, Denmark office:
Italy Switzerland
E-mail: opencadd@sanet.com.br Tel: +45 (0)39-66 56 50
Teoresi s.r.l. Scientific Computers SC AG
Authorized Reseller Fax: +45 (0)39-66 56 20
Tel: +39 (0)11-240-80-00 Tel: +41 (0)31-954 20 20
E-mail: info@comsol.dk
Czech Republic, Slovakia, Russia, Fax: +39 (0)11-240-80-24 Fax: +41 (0)31-954 20 22
Web: www.comsol.dk
Ukraine, Belarus, Moldavia E-mail: info@teoresi.it E-mail: info@scientific.ch
Web: www.teoresi.it Helsinki, Finland office: Authorized Distributor
Humusoft s.r.o.
Authorized Distributor Tel: +358 (0)9-455-00-55
Tel: +420-2-68-44-174
Fax: +358 (0)9-455-00-51 Taiwan
Fax: +420-2-68-44-174
Japan E-mail: info@comsol.fi Scientific Formosa, Inc.
E-mail: byron@humusoft.cz
Cybernet Systems Co., Ltd. Web: www.comsol.fi Tel: +886 (0)2-2505-0525
Web: www.humusoft.cz
Tel: +81 (0)3-5978-5410 Trondheim, Norway office: Fax: +886 (0)2-2502-4478
Authorized Reseller
Fax: +81 (0)3-5978-5440 Tel: +47 (0)735-09-220 E-mail: info@pent1.greenhills.com.tw
France E-mail: infomatlab@cybernet.co.jp Fax: +47 (0)735-09-221 Authorized Distributor
Scientific Software Group Web: www.cybernet.co.jp E-mail: info@comsol.no
Authorized Distributor Web: www.comsol.no U.K., Ireland
Tel: +33-01-41-14-67-14
Authorized Distributors Cambridge Control Ltd.
Fax: +33-01-41-14-67-15
Korea Tel: +44 (0)1223-423-200
E-mail: info@ssg.fr
Kimhua Technologies, Inc. Poland Fax: +44 (0)1223-423-289
Web: www.ssg.fr
Tel: +82 (0)2-556-1257 ControlSoft E-mail: info@camcontrol.co.uk
Authorized Distributor
Fax: +82 (0)2-556-4020 Tel: +48 (0)12-6-17-33-48 Web: www.camcontrol.co.uk
Germany, Austria E-mail: info@soft.kimhua.co.kr Fax: +48 (0)12-6-33-27-12 Hove, England office:
Scientific Computers GmbH Web: kimhua.co.kr E-mail: info@controlsoft.krakow.pl Tel: +44 (0)1273-722-838
Tel: +49 (0)241-470-750 Authorized Distributor Web: www.controlsoft.krakow.pl Fax: +44 (0)1273-720-550
Fax: +49 (0)241-449-83 Authorized Distributor E-mail: pauline_fox@camcontrol.co.uk
The Netherlands, Belgium, Luxembourg
E-mail: matlab.info@scientific.de Authorized Distributors
Scientific Software Benelux B. V. Singapore, Malaysia, Thailand,
Web: www.scientific.de
Tel: +31-(0)182-53-76-44 The Philippines, Indonesia, Brunei
Unterföhring (Munich) office: Fax: +31-(0)182-57-0380 TechSource Systems Pte Ltd.
For inquiries outside the U.S.
Tel: +49 (0)89-995-901-0 E-mail: info@ssb.nl and Canada, please contact your
Tel: +65-842-4222
Fax: +49 (0)89-995-901-11 Web: www.ssb.nl
Authorized Distributors
Fax: +65-842-5122 local distributor. If your country
Authorized Distributor E-mail: info@techsource.com.sg is not listed, please contact
Authorized Distributor
The MathWorks directly.
Contacting The MathWorks

Please send e-mail to info@mathworks.com, fax 508-647-7101, call 508-647-7000, or access our
Web site at www.mathworks.com. For address changes, please send e-mail to access@mathworks.com.
Useful e-mail addresses:
info@mathworks.com Sales, pricing, and general information
support@mathworks.com Technical support
access@mathworks.com MATLAB Access program information
suggest@mathworks.com Product enhancement suggestions
resumes@mathworks.com Resume submission for job opportunities
news-notes@mathworks.com MATLAB News & Notes Editor
updates@mathworks.com Microcomputer updates and subscriptions
service@mathworks.com Customer Service: order status, license renewals, passcodes
finance@mathworks.com Financial products information
doc@mathworks.com Documentation error reports
bugs@mathworks.com Bug reports
connections@mathworks.com Third-party products and services
Internet services:
www.mathworks.com The MathWorks home page
ftp.mathworks.com FTP server
Other resources:
comp.soft-sys.matlab Usenet newsgroup 24 Prime Park Way
Natick, MA 01760-1500 USA
© 1998 by The MathWorks, Inc. All rights reserved. MATLAB, Simulink, Handle Graphics, and Real-Time Workshop are registered trademarks, and Stateflow and
Target Language Compiler are trademarks of The MathWorks, Inc. Other product or brand names are trademarks or registered trademarks of their respective holders.
9521v00 2/98

Using Simulink in Engine Model PDF

Uploaded by

Copyright:

Available Formats

Using Simulink in Engine Model PDF

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Using Simulink in Engine Model PDF

Uploaded by

Copyright:

Available Formats

S IMULINK -S TATEFLOW T ECHNICAL E XAMPLES

Using Simulink and Stateflow

System Models in Simulink................................................................................................ 7

System Models in Simulink with Stateflow Enhancements .......................................... 49

USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 3

hydcyl4.mdl — Four-cylinder model

hydrod.mdl — Two-cylinder model with load constraints

4 S IMULINK -S TATEFLOW T ECHNICAL EXAMPLES

VII. Automatic Transmission Control

VIII. Electrohydraulic Servo Control

IX. Modeling Stick-Slip Friction

Andy Grace agrace@mathworks.com

Paul Barnard pbarnard@mathworks.com

Larry Michaels lmichaels@mathworks.com

Bill Aldrich baldrich@mathworks.com

USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 5

Attention: Robert W. Weeks

Attention: Sham Ahmed

6 S IMULINK -S TATEFLOW T ECHNICAL EXAMPLES

USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 7

8 S IMULINK -S TATEFLOW T ECHNICAL EXAMPLES

R = specific gas constant

Intake Mass Flow Rate

USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 9

N = engine speed (rad/s)

Torque Generation and Acceleration

Torque eng = −181.3 + 379.36 m a + 21.91( A / F ) − 0.85( A / F ) 2 + 0.26σ − 0.0028σ 2

ma = mass of air in cylinder for combustion (g)

JN˙ = Torqueeng − Torqueload Equation 1.5

J = Engine rotational moment of inertia (kg-m2)

10 S IMULINK -S TATEFLOW T ECHNICAL EXAMPLES

Engine Timing Model in Simulink 2

mass(k) Air Charge

Figure 1.1: The top level of the Simulink engine model

USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 11

Throttle Manifold Dynamics

Throttle Ang. Limit to Positive

1.0 Atmospheric Pressure, Pa (bar)

mdot Input (g/s) mdot to Cylinder (g/s) 1

Throttle Flow vs. Valve Angle and Pressure

Intake Manifold Vacuum

Figure 1.2: The Throttle and Intake Manifold Subsystems

12 S IMULINK -S TATEFLOW T ECHNICAL EXAMPLES

USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 13

Figure 1.3: The simulated engine speed

Modeling — SPEED C ONTROL

14 S IMULINK -S TATEFLOW T ECHNICAL EXAMPLES

mass(k) Air Charge

Figure 1.4: A discrete-time PI controller is added to the engine model

USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 15

Figure 1.5: Speed Controller Subsystem

Closed–Loop Engine Speed Control

Figure 1.6: Typical simulation results

16 S IMULINK -S TATEFLOW T ECHNICAL EXAMPLES

USING S IMULINK AND S TATEFLOW IN A UTOMOTIVE A PPLICATIONS 17

Two DOF Spring/Damper Model

Figure 4.3: The Spring/Damper suspension subsystem

Lf = 0.9; % front hub displacement from body CG

[t,x] = sim('suspn2',10); % run a time response