Volume 4 No.1, Winter 2000

ISSN# 1523-9926


Complex Systems Simulation Using


Mikhail M. Polonskii
Universidade de Passo Fundo  - RS


A procedure for complex systems simulation using MATLAB/SIMULINK software is presented. As an example,  an internal combustion engine was chosen. The procedure is explained in detail to help for university students and industry specialists in mastering  this software. The models (MATLAB .m and .mdl-files) developed are available for download.


The use of graphical dynamic system simulation software is becoming more popular as  engineers strive to reduce the time to develop new control systems [ 1 ]. Dynamic system simulation software is an important tool for developing advanced reliable and high quality products and systems. Thus, the skill of complex systems simulation is the inherent goal of any engineer curriculum. In this article we show how to convert a set of differential non-linear time variant equations into easy to use MATLAB/SIMULINK graphical model. MATLAB/SIMULINK software has been chosen because many people use this software in both academia and industry. As the object of simulation we have selected an internal combustion  engine.

The use of model based control methods designed to meet future emission and diagnostic regulations has also increased the need for validated engine models. A previously validated, nonlinear, mean torque predictive engine model is converted to MATLAB / SIMULINK to illustrate the benefits of a graphical simulation environment. The model is used to simulate a sequencial-fuel-injected, spark ignition engine and includes air and fuel  dynamics in the intake manifold as well as the process delays inherent in a four stroke cycle engine. The engine model can be used in five ways [ 1]:

Although developed  for a specific engine, the model is generic enough to be used for a wide range of spark-ignition engines. Modular programming techniques reduce model complexity by dividing the engine  into hierarchical subsystems. Further discussion supposes some basic knowledge of MATLAB/SIMULINK of the reader.


The engine model for a four stroke spark ignition engine was developed by Moskwa and Hendrick [ 2 ]. It was subsequently refined and used for different studies. Here we use the model as it has been presented by Godbole [ 3, 4 ]. This is continuous time model with three state variables given by:
    - the mass of air in intake manifold;
    - the engine speed;
    - mass flow rate of fuel entering combustion chamber.
The exhaust gas recirculation (ERG) system has not been included in this model. The following state equation is obtained by applying law of conservation of mass to the air flow in the intake manifold:

                      ( 1 )

where  = mass of air in the intake manifold;    = mass rate of air in the intake manifold;
= mass rate of air entering the intake manifold; = mass rate of air leaving the intake manifold and entering the combustion chamber.
The mass of air entering the intake manifold is modeled as follows:

           ( 2 )

where  = the maximum flow rate corresponding to full open throttle,  = Normalized Throttle characteristic,
= Normalized Pressure Influence function.
Usually  is determined by an experiment and a data table is used in the simulation. Here we use a curve fit as in [ 2 ]:

              ( 3 )

where  = the throttle angle.
The Normalized Pressure Influence function  is calculated as a function of manifold to atmosphere pressure ratio:

               ( 4 )

where  = intake manifold pressure,  = atmospheric pressure.
The intake manifold pressure  and the mass of air  entering the intake manifold are modeled by ideal gas law:

                                ( 5 )

where  = constant value,  = constant value,  = constant value,  = intake manifold volume

The flow of air from intake manifold to combustion chamber is given by

                      ( 6 )

where  = engine angular velocity,  = volumetric efficiency,  is given by

                                         ( 7 )

The volumetric efficiency  represents the effectiveness of engine's induction process. It is a complex function of engine geometry and many engine parameters. It is normally modeled by experimental data. Here  the empirical formula given in [ 2 ]
is used:

                   ( 8 )

The dynamics of the fuel injection process is modeled as follows:

                      ( 9 )

where  = fuel rate entering the combustion chamber,  = command fuel rate,  = effective fueling time
constant, which is modeled as

            ( 10 )

where  = desired air/fiel ratio (tuning parameter),  = constant value.
The third state equation is obtained by applying Newton's second law to the rotational dynamics of the engine:

             ( 11 )

where  = effective inertia of the engine,  = engine indicated torque,  = engine friction torque,  = accessories
In a four stroke engine, the process of torque production is discrete depending on the engine speed. The model we used here is continuous, thus two delays were introduced. The indicated torque is modeled by:

      ( 12 )

where,  = intake to torque production delay,  = spark to torque production delay,  = normalized air fuel
ratio influence function,  = normalized spark influence function,  = the maximum torque production capacity of an engine for  = 1  and  = 1.  Here the following curve fit relation from [ 2 ] is used:

               ( 13 )

where  = actual air/fuel ratio of the mixture in the combustion chamber. With some simplifying assumption, it will be given by:

                       ( 14 )

Here  the following curve fit relation from [ 2 ]is used:

           ( 15 )

where  = tuning parameter representing the spark advance from Top Dead Center (TDC),  = the minimum spark advance from TDC for best torque. The engine friction torque is curve fitted from experimental data [ 2 ] as

          ( 16 )

Thus the engine model has tree states  (    ,  ), two control parameters (  ) and one control variable  (  ).  Table 1 summarizes the parameter values used for simulation, which had been obtained for a 6-cilinders, 2.7 liter car engine.                                                                                    

0.1843  kg/s

0.0038  m^3

0.0027 m^3

0.1454  kg m^2

498636  N m/ (kg/s)

5.48/     s

1.3/     s

300  deg  K

28.84  g/mole

8314.3  J/ mole deg K

Table 1. Simulation Parameter Values


The complete model of the engine can be divided into air,  fuel and  torque subsystems. This division is logical because the model is based on the three state variables ( ). The air subsystem involves the equations from Equation 1 to Equation 8. The fuel subsystem is determined by  Equation 9 and Equation 10, and the torque subsystem is defined by the
equations from Equation 11 to Equation 16.
To create the MATLAB/SIMULINK model of the engine, let us begin with the air subsystem. The model should be hierarchical because this simplifies their use and modification. Thus, we'll build the model through forming the subsystems while the subsystems will be created by grouping simple blocks. First, we'll transform Equations 6, 7 and 8 into MATLAB/SIMULINK model. Analyzing these equations, it's easy to deduce that the inputs are  (  , )   and the output is ( ). Three inputs for this subsystem are the model parameters (  ) and these parameters will be introduced as constants at the top level of hierarchy of the model.  Equations 6, 7 and 8 determine the air_3 subsystem of the air subsystem as shown in Figure 1.

Figure 1. The simulation diagram of the  air_3 subsystem.

After creating this subsystem, we group all the blocks and entitle this subsystem as air_3.  The input and output ports are labeled as the corresponding variables. The next step is to convert  Equation 5 into MATLAB/SIMULINK model. It is important to realize that  depends on  , thus we'll get the air_2 subsystem as follows:

Figure 2. The simulation diagram of the air_2 subsystem.

The last part   of the air subsystem model is air_1  and this involves  Equations 2, 3 and 4. The MATLAB/SIMULINK corresponding model  is present in Fig. 3.

Figure 3. The air_1 subsystem's MATLAB/SIMULINK model.

Now we can create the complete air subsystem model. This model joints the air_1, air_2, air_3 subsystems and involves one integrator and one sum block as appears in Fig. 4.

Figure 4. The complete air subsystem  model.

Now we create the torque subsystem model. The model involves  Equations 11, 12, 13, 14, 15 and 16. Obviously, this subsystem is complex and should be divided  into some subsystems. The first subsystem is torque_1. This subsystem is based on the Equations 13, 14, 15 and calculates  for Equation 12.  Figure 5 shows this model.

Figure 5. The torque_1 subsystem model.

This model involves four variable transport delay blocks  according to Equation 12. The point is to adjust the initial input value for each of these blocks. The initial input for the Variable Transport Delay block (Figure 5) is specified as 100 which is the initial engine velocity value. The initial input for the Variable Transport Delay 1 block was adjusted to 0.014 because this value had been obtained as the steady-state value for the velocity equals to 100 1/s. The initial input for the Variable Transport Delay 2 block was adjusted to 0.5 which is a reasonable value for cosine function. The initial input for the Variable Transport Delay 3 block was adjusted to 0.5  too. These initial  specifications determine the initial part of transient of the engine velocity.
The presence of the Constant, the Relational Operator and the Switch blocks (Figure 5) should be commented. If we send the output of the  Variable Transport Delay 2 block directly to the corresponding input of the Mux block, the  transient  of the engine velocity at the beginning of simulation will tend to fall down. To discover the reason of such strange model dynamics we can connect one Auto-Scale Graph ( for MATLAB v.4.2 ) to output of  the Variable Transport Delay 2 block. Then we can see that the output of the Variable Transport Delay 2 block is oscillatory around the zero value. This output is used as the factor for indicated torque calculation (Equation 12). Thus, using  the Constant, the Relational Operator and the Switch blocks (Figure 5) we can clip the output of the Variable Transport Delay 2 block to zero if this output becomes negative.
The torque_2  subsystem was created to simulate  Equations 11,12 and 16.  Figure 6 shows the model.

Figure 6. The torque_2 subsystem's model.

 It is necessary to define the initial value of the Integrator (Figure 6) equal to 100, which is the initial value of the engine speed. If the initial value of the Integrator is zero, the MATLAB zero-devision error  is generated.
The complete torque subsystem is shown in Figure 7

Figure 7. The torque subsystem  model.

The last subsystem is the fuel subsystem and this involves  Equations 9 and 10. First, we create the   FUEL_1 subsystem model based on  Equation 9 as shown in Figure 8.

Figure 8. The MATLAB/SIMULINK FUEL_1 subsystem model.

This realization is better than one with numeric differentiation because of the absence of algebraic loop. The complete model of the fuel subsystem is done in Figure 9.

Figure 9. The model of the fuel subsystem.

Figure 10 shows the complete engine model with all the links between the subsystems.

Figure 10. The whole model of the engine.

We can group these three subsystems into one and thus get the final version of the model. This version is shown in Figure 11.

Figure 11. The final version of the engine MATLAB/SIMULINK model.

 All of the simulation parameters  are present at this upper level of the model, which simplifies the model modification. The values of the parameters    and    have been adjusted for the best engine torque production. Simulation of the variations of these parameters can be useful to fuel injection and ignition systems projects.
This model is available for download as  eng.m (for MATLAB v. 4.2) or as eng.mdl (for MATLAB v.5).


Different types of simulation can be done using this MATLAB/SIMULINK engine model. Here we simulate the open-loop dynamics of the engine. The first simulation is executed  for the engine when the external torque is zero, the throttle angle is constant  = 1 rad (Figure 11) and the engine velocity transient response is obtained. The result is presented in Figure 12.

Figure 12. The engine velocity output (rad/s) in response to a unit throttle position input.

The velocity transient response shows some time delay at the beginning of the process and the steady state velocity is about 620 rad/s. The initial time delay reflects the presence of the time delay factors in  Equation 12.
Another simulation is done for  the throttle position  = 0.6 rad  when the external torque = 10 sin t.  The engine velocity transient response is shown in Figure 13.

Figure 13. The engine velocity output,  = 0.6 rad, external torque10 sin t.

This last simulation shows that the steady state velocity is about  370 rad/s and there are some velocity oscillations according to the extern torque variations. Thus, to reduce the output sensitivity to extern torque variations, it is necessary to implement some type of feedback control scheme.
Now we'll see how the engine model can be used for closed-loop simulations. Let us investigate the dynamics of a idle speed
relay controller. The model is shown in Figure 14.

Figure 14. Idle engine speed control system

Here some of the engine parameters were hidden inside the ENGINE subsystem. The controller controls spark advance angle, which is being switched between 0.2 and 1, i.e. between 20% and 100% of the best torque spark advance. During idle mode the engine works without load and their response is very fast. This is the reason to control the spark advance, instead of the throttle angle, because  is more than 4 times greater than  (Table 1). The throttle angle  = 0.32  rad, it has been tuned by simulations and is contant. The REFERENCE  SPEED  signal is shown in Figure 15.

Figure 15. The REFERENCE  SPEED signal.

Initially the REFERENCE SPEED is 150 rad/s and is greater  than the idle speed of the engine, which is 100 rad/s. If we try to simulate the engine response, when the REFERENCE SPEED is 100 rad/s initially, the engine stops due to engine friction torque and insufficient controller output.
The External Torque signal is shown in Figure 16. This signal has been chosen to verify the system response under rapid external torque variations.

Figure 16. The External Torque signal.

The SPARK ADVANCE signal obtained by simulating the model (Figure 14) is shown in Figure 17.

Figure 17. The SPARK  ADVANCE  signal.

It can be seen that the frequency of the relay switching is variable and depends on tne external torque value.
The ENGINE SPEED response is shown in Figure 18.

Figure 18. The ENGINE SPEED response.

We can see, that the engine response is quite fast and precise. The system is not sensitive to external torque variations. The procedure of the simulation has been as follows:
    1. The open-loop dynamics of the engine has been simulated to define the convenient constant  value;
    2. The ralay controller has been added and switching limits has been chosen;
    3. Variable torque has beed added and final adjustment of the ralay  has been done.

The model ( Figure 14 ) is available for download as   idle.m  (for MATLAB v. 4.2).


This paper has presented the approach to complex systems simulation using MATLAB/SIMULINK. As an example, the equations of the four-stroke, 6-cilinder, 2.7 liter internal combustion engine has been chosen. The "down-up" method of hierarchical model creation has been explained step-by-step. Some simulations have been executed to demonstrate the open-loop and closed-loop dynamics of the engine. The model can be used for closed-loop cruise velocity control and injection/ignition systems simulations, too. It is possible to integrate this model into the complete model of car dynamics [ 4 ]. The results presented here can be adopted to the electronics, mechatronics and mechanics curricula to help in development of the skill of complex system simulation.


[ 1 ] R.W. Weeks and J.J. Moskwa, "Automotive Engine Modeling for Real-time Control Using MATLAB/SIMULINK", SAE paper No. 950417, 1995.

[ 2 ] J.J. Moskwa and J.K. Hedrick, "Modeling and validation of automotive engines for control algorithm development", ASME Journal of Mechanisms, Transmissions and Automation in Design, vol. 114, no. 2, pp. 278-285, 1992.

[ 3 ] D. Godbole and S. Karahan, "Automotive Powertrain Modeling, Simulation and Control using Integrated System's CASE Tools". Integrated Systems, Inc. Santa Clara, CA.

[ 4 ] D. Godbole, "Automotive Demos for Xmath. Automotive Suspension Design and Powertrain Speed Control". Intelligent Machines and Robotics Laboratory. University of California, Berkeley.   godbole@robotics.eecs.berkeley.edu

* MATLAB and SIMULINK are trademarks of The MathWorks, Inc.


Return to This Issues Home Page