High-fidelity simulations of dynamic embedded systems can be invaluable. This follow-up to “Modeling Dynamic Systems” (August 2000) presents some techniques and algorithms you might find useful.
By Jim Ledin
In a previous article, I discussed approaches for developing mathematical models of dynamic system components. A dynamic system has behavior that evolves over time, often in a manner that can be represented by differential equations. Our interest in dynamic systems is in those that contain embedded computing hardware and software, which we will refer to as dynamic embedded systems. Some examples of dynamic embedded systems are automobile engine control and cruise control systems, assembly line robots, satellite attitude control systems, air-to-air guided missiles, and Mars rovers.
Although mathematical models of dynamic system components can be used in several different ways to analyze and predict system behavior, our interest here is in assembling a set of models and executing them over time as a simulation.
The use of simulation in the development of dynamic embedded systems goes back to the earliest days of embedded computing, when digital processors were first designed into military weapons. Simulation is useful at many stages of the product development process, from initial concept exploration through the testing of maintenance software releases. Simulations can easily include complex real-world non-linearities that may be ignored or simplified by other approaches such as linear systems analysis.
One particularly useful type of simulation is hardware-in-the-loop (HIL) simulation, which involves connecting a real-time simulation to an embedded computing device. HIL simulation testing enables thorough, repeatable testing of significant parts of an embedded system across their operational envelopes in a laboratory environment. This approach typically leads to a more thoroughly tested product at significantly lower cost than with traditional test approaches.
To develop a simulation of a dynamic system, we must first construct mathematical models of suitable fidelity to describe the system components and their interactions with the environment. We then assemble the system simulation from these models and other elements such as integration algorithms. The models used in a simulation must be mutually compatible so that signals flow between the models with the proper units and timing.
A complete simulation must include the correct initialization of state variables and other parameters during simulation start-up, as well as evaluation of the model dynamic equations and integration of state variables during simulation execution. In addition, the simulation must have a user interface to allow the user to specify initial conditions and model parameters and to control the operation of the simulation.
In simulating a continuous dynamic system on a digital computer, we estimate the solution of the differential equations that describe the system using a numerical integration algorithm. Many different kinds of numerical integration algorithms are available. We must select one that is appropriate for our simulation.
Numerical integration algorithms
We will assume that the behavior of a continuous dynamic system is defined by a set of coupled first order differential equations. These equations are frequently nonlinear, which makes their analytic solution difficult or impossible. As a result, we must resort to numerical algorithms to estimate the solution of these differential equations.
We will require that all differential equations that we work with be of first order since our integration algorithms will only work with first order equations. We can transform an arbitrary differential equation containing second (or higher) derivatives into an equivalent set of first order equations.1
We will discuss numerical integration algorithms in terms of the differential equation:
where f(x,u(t)) is a possibly nonlinear function and u(t) is an arbitrary input function.
One of the simplest algorithms for performing numerical integration arises from the definition of the derivative of a continuous function, which is shown in Equation 1:
If we use a small but finite value of h, we can approximate the derivative as shown in Equation 2:
In terms of sampling times, with the integer subscript n indicating t = nh and with a hat over a variable indicating an approximate solution, we now have Equation 3:
This is the forward Euler integration algorithm. It is a first order integration method, because it contains the function f evaluated at one point in time. It is also an explicit method, which means that it does not use states or inputs from future values of n. This algorithm is shown graphically in Figure 1, where the cross-hatched area is the approximation of the integral of f over the interval between t = nh and t = (n+1)h.
Alternatively, we can use a backward difference to approximate the derivative, as shown in Equation 4:
This results in the difference equation shown in Equation 5, which is known as the backward Euler integration algorithm:
This algorithm is shown graphically in Figure 2.
The backward Euler integration algorithm is a first order implicit method. It is implicit because the state is a function of itself. Implicit algorithms require additional computation to solve for , which commonly takes the form of an iterative technique such as the Newton-Raphson algorithm. Implicit algorithms have advantages in terms of accuracy and numerical stability over explicit methods. However, the drawbacks of implicit algorithms are the additional computation requirements and their inappropriateness for real-time applications.
Implicit methods are not suitable for real-time use for two reasons. First, the execution time required for the iterative solution of is unpredictable. Second, an input value from a future time is required in the form of un+1, which is not available at time step n when the integration step must be performed.
A wide variety of other types of numerical integration algorithms are available, many of which possess unique attributes that are valuable in specific applications. For example, the Runge-Kutta algorithm family is particularly good at simulating systems containing discontinuities in the function f. To achieve acceptable accuracy, integration algorithms of second through fourth order are commonly used instead of the first order algorithms we discussed above. This is because higher-order methods tend to produce more accurate results than the lower order methods. This leads us into a discussion of the sources of error in numerical integration algorithms and ways to minimize that error.
Since the difference equations that implement these integration algorithms necessarily produce approximations of the exact solution of the differential equations, each step of an integration algorithm will introduce some error into the estimated solution. The error that occurs in each integration step arises from two sources: roundoff error and truncation error.
Roundoff error is the result of a computation that uses finite-precision rational numbers (such as the floating-point numbers used by computers) to approximate infinite-precision real numbers. The amount of roundoff error in a particular computation depends on the number and type of arithmetic operations performed and the precision of the floating-point representation used. Given a specific computer and sequence of arithmetic operations, the only way we can reduce roundoff error is by increasing the precision of the number representation used in the computation. We can do this, for example, by switching from single precision floating-point to double precision. Single precision floating-point numbers typically have about seven decimal digits of precision, while double precision floating-point numbers have about 15. The drawbacks in using double precision rather than single precision are that the double precision variables take up twice as much space in memory and more processor time may be required to perform the computations.
For simulations of practical systems, the use of single precision floating-point numbers (rather than double precision) may introduce significant roundoff error into the solution of the dynamic equations. The severity of this problem depends on the form of the dynamic equations as well as the integration algorithm selected, the integration time step size, and the length of the simulation run. As a rule, when performing numerical integration, it is preferable to use double precision in the integration algorithms whenever possible. The additional memory and processing time required by the use of double precision are rarely significant costs in comparison to the benefits of increased simulation accuracy.
The second source of error, truncation error, is the result of an integration step wit no roundoff error. In other words, truncation error is the error if the integration algorithm step were executed on a computer with infinite numberical precision.
Each integration algorithm has its own truncation error characteristics. If we assume that the function f is continuous over the integration step and we are using a sufficiently small step time h, we can estimate the truncation error for a single integration step using Equation 6. Here, the constant a depends on both the local characteristics of the function f and the integration algorithm used. k is the order of the integration algorithm and h is the integration step size:
From Equation 6 we see that reducing the integration step size can reduce the truncation error to an arbitrarily small level. However, doing this will increase the execution time for the simulation while increasing the roundoff error in the solution. Using a higher order integration algorithm will also reduce the truncation error, but this approach has drawbacks as well. A higher order integration algorithm will require additional execution time, but in most cases the primary drawback of the higher-order algorithms is that the possibility of numerical instability increases.
Integration algorithm stability
Numerical instability occurs when roundoff and truncation errors, which may be very small during initial integration steps, grow without bound over a number of steps. Integration algorithm instability occurs when the time step size is too large, which causes the difference equations used to estimate the solution of the dynamic equations to become unstable.
Each combination of dynamic system and integration algorithm has its own stability characteristics. In our discussion here, we will assume that it is possible to linearize a simulated nonlinear system about an equilibrium point. With this assumption, it will always be true that given a small enough time step size, integration algorithm stability is achievable as long as the dynamic system being simulated is stable near the equilibrium point.
We will discuss numerical stability in terms of the linear first order dynamic system shown in Equation 7, with the exact solution as given in Equation 8:
We can examine the effect of various integration time step sizes on the results of the numerical integration. Figure 3 shows the solutions of Equation 7 using the forward Euler method of Equation 3 with six different values of the time step h. Figure 3a, with a value of 0.01 for h, has good accuracy. Figure 3b has h = 0.5 and there is significant error, though the trend of the solution is still in reasonable agreement with the exact solution. Figure 3c has a value of 1 for h and the solution goes to zero after one step and stays there, which is clearly inaccurate, though stable. Figure 3d has h = 1.9 and the solution is a damped, stable oscillation. Figure 3e has h = 2 and the solution is an undamped oscillation with constant amplitude. Figure 3f has a value of 2.1 for h and the solution is an oscillation with growing magnitude, which represents an unstable system. Larger values of h will produce solutions with magnitudes that grow at faster rates.
For this simulation using the forward Euler integration method, any value of the time step size h greater than 2 produces a numerically unstable solution. As we can see from Figure 3, by the time that h has increased to the point of instability, any concept of accuracy in the solution has long been lost. This will happen in situations where the time constants of the dynamic system elements have values that do not vary from each other by more than a few orders of magnitude. The other case consists of systems with time constants that vary over several orders of magnitude. These are called stiff systems.
A system may include both fast-moving and slow-moving dynamics, which are represented by the time constants of the system model. We characterize the stiffness of a system as the ratio of its longest time constant to its shortest time constant. We consider a system to be stiff if this ratio is roughly three orders of magnitude or higher. When simulating stiff systems, we must use care in the choice of integration method to ensure that the solution remains numerically stable and produces accurate results while executing with reasonable efficiency.
Consider the open-loop dynamic system represented by the Simulink model shown in Figure 4. This system consists of a first order actuator driving a first order plant. In this example, the actuator has a time constant of 1ms and the plant has a time constant of 1s. The system is driven by a unit step input.
If we simulate the system of Figure 4 using the forward Euler integration method and a step time of 1.9ms for a period of 10s, we find that the peak magnitude of the integration error is about 0.001. If we increase the step time slightly to 2.01ms the simulation becomes numerically unstable and the truncation error grows without bound. The point of this example is that even in this case of mild stiffness, there is only a little degradation of simulation accuracy before instability occurs.
Many real world dynamic systems have several orders of magnitude more stiffness than the example of Figure 3. To simulate these very stiff systems with an arbitrarily selected integration algorithm, it may be necessary to take extremely small time steps to accommodate the shortest system time constant. This is inefficient because the system components with longer time constants do not require the same small integration steps needed to simulate the faster components. In fact, if a very small time step is used for the slow-moving components, the accuracy of the overall simulation may suffer due to roundoff error that accumulates over the large number of integration steps required.
Two approaches are commonly used in the simulation of stiff systems. The first is to select an integration algorithm that works well in the presence of stiffness. For non-real-time simulation applications, many integration algorithms are available that perform well with stiff systems. The implicit algorithms introduced above, for example, have good performance in the presence of stiffness. The Gear family of implicit algorithms is intended primarily for the simulation of stiff systems.3 By using integration algorithms with good performance in the presence of stiffness, it is possible to maintain simulation accuracy and stability while taking much larger time steps than would be possible with other algorithm types such as Runge-Kutta.
A second approach to the simulation of stiff systems involves the use of multiframing. Multiframe simulations break the simulation into multiple parts that execute at different frame rates. In the example of Figure 4, we might run the actuator simulation with a 0.1ms step time and the plant simulation with a 0.1s step time.
The multiframe method adds complexity because we must extrapolate or interpolate any values from the slower frame rate subsystem that are used as inputs by the fast frame rate subsystem. There may also be some difficulty in performing I/O at fixed intervals if the various frame rates are running on a single processor.
Of these two approaches to simulating stiff systems, only the multiframing method is applicable for real-time applications. This is due to the unsuitability of implicit integration algorithms in situations where real time I/O is required. In non-real-time simulations, it is often simplest to employ an integration algorithm with good performance in the presence of stiffness, which avoids the difficulties associated with multiframing.
Combined discrete-continuous systems
A dynamic system modeled with both difference equations and differential equations is called a combined discrete-continuous system, or simply a combined system. Combined systems present special challenges for simulation because the integration time step required for the simulation of the continuous system components often varies from the time step of the discrete subsystem. In addition, the inputs from the discrete components to the continuous components may appear as discontinuous signals, which limits the selection of integration algorithms.
Figure 5 shows a common situation in which a continuous plant is controlled by a digital controller. The controller operates at a fixed update rate and receives an input voltage from the plant via an analog-to-digital converter that samples and quantizes the signal. Outputs from the digital controller are sent to a digital-to-analog converter, which converts the numerical value to a voltage and holds it at a constant value until the next update.
In the most general situation, where the optimal integration time step for the continuous plant simulation differs from the controller update interval, the simulation must use the multiframe technique. For an accurate simulation, the plant model will normally require a smaller integration time step than the controller update interval. When using a fixed-step integration algorithm, the simplest approach is to have the integration step size be an integer submultiple of the controller update interval.
When the inputs from the slower rate system to the faster rate system are digital-to-analog converter outputs, no extrapolation is required because the signals remain constant between updates. This simplifies the interface between the two frame rates somewhat.
Combined systems are an important application area for dynamic system simulation. Simulation is often used in the design, development, and test of complex dynamic systems in which an analog plant is controlled by an embedded digital controller.
We will now discuss the design and implementation of a complete simulation of a dynamic system. The system we choose to simulate is a rotating turntable driven by a DC motor. A digital controller accepts the commanded angular position for the turntable as an input and generates a voltage to drive the DC motor. The motor moves the turntable to the commanded position. This is a combined discrete-continuous system which we will simulate using the multiframing technique in Simulink.
We model the controller as a digital system that operates at an update interval of 5ms. The input to the controller is the turntable position as measured by a quadrature position encoder. This device outputs a digital angular position reading with a resolution of 4,096 steps per turntable revolution. The controller output signal is converted to a voltage via an 8-bit digital-to-analog converter (DAC), which drives the DC motor through a power amplifier that has a voltage gain of 1. We wish to use as inexpensive a processor as possible, so we will model the controller calculations as taking 4ms to complete during each update cycle.
When we model a dynamic system, we must decide how much detail to include in the representation of each component. The amount of detail required will depend on the intended uses for the simulation. In this case, the primary goal of the simulation is to enable the design and implementation of a digital controller for the turntable system. Consequently, we will attempt to model the motor and turntable dynamics with good accuracy and we can ignore the dynamic behavior of the power amplifier. We will want to model the behavior of the digital control algorithm as accurately as we can.
Figure 6 shows a Simulink block diagram of this system. The motor dynamics are modeled as a second order transfer function. The input to this transfer function is the motor drive voltage and the output is the motor velocity. The motor velocity is integrated to determine the motor position. The position encoder outputs a quantized measurement of the motor position to the digital controller. The inputs to the digital controller are the commanded turntable position and its current position as measured by the encoder. The controller outputs an analog voltage that drives the power amplifier, which in turn drives the motor. The boxes labeled “DAC Output” and “Position (radians)” represent signals that are saved during simulation execution for later plotting.
Looking inside the box in Figure 6 labeled “Digital Controller,” we see the implementation of the control algorithm, which appears in Figure 7. Since this algorithm is to be executed at intervals of 5ms, Simulink requires that we add the blocks labeled “Zero-Order Hold1,” “Zero-Order Hold2,” and “Unit Delay.” Each of these blocks has a parameter associated with it that has been set to indicate 5ms updates. The remainder of the algorithm consists of a proportional-derivative controller that uses the discrete-time transfer function labeled “Derivative Estimate” to produce an estimate of the turntable angular rate. The output of the control algorithm is fed to an 8-bit DAC that limits and quantizes the signal. The “4ms Transport Delay” block delays the output of the algorithm by 4ms, which simulates a computational delay.
In addition to implementing the models for the simulation, we must also specify the integration algorithm to use and (if the algorithm is of the fixed-step type) the integration time step. After some experimentation, I determined that a 1ms time step size combined with Simulink's fourth-order Runge-Kutta integration algorithm provide satisfactory results. It is easy to confirm this by performing a run with a step time of 0.1ms and verifying that the results do not change significantly.
Figure 8 shows the step response of this system. The initial state of the system is with the turntable stopped at the zero position. At time zero, a command is applied to move the turntable to the 1 radian position. We can see that by about 5s, the turntable has essentially reached the desired location. This plot gives little indication that the nonlinearities we have modeled have any effect on the system.
We can examine the output of the digital controller as shown in Figure 9. Here, we can clearly see the effects of the quantization of the position encoder signal. We can also see that the DAC output signal is saturated at the maximum value for about the first 1.7s of the simulation and the quantization of the DAC output signal voltage is visible as well.
This simulation of the turntable and controller provides us with a great deal of information about the dynamic behavior of the overall system. The details of nonlinear system behavior are crucial when developing designs that must work in the real world. The use of simulation allows us to determine the behavior of the system under a wide variety of conditions by changing model parameters and by adding additional effects and complexity to the component models. We can also easily perform trade studies to determine the behavior of the system in response to design modifications, such as a different controller update rate, different number of DAC bits, or different encoder resolution.
We have discussed several of the techniques and issues involved in the development and use of simulations of dynamic embedded systems. We also examined a simulation of a dynamic system that included a number of realistic, nonlinear effects. We saw how the use of simulation enables us to analyze the detailed behavior of embedded control algorithms before implementing them in an actual system.
Simulation is an important and commonly used technique in the development and testing of complex embedded systems. We can use simulation beginning with product concept definition through development engineering and on to the testing of maintenance software releases and hardware upgrades. Appropriate use of simulation in a development project can reduce the overall time and cost required for the project while producing a more thoroughly tested, higher quality product.
Jim Ledin , PE, is a consulting engineer in Camarillo, CA. He is the author of Simulation Engineering (to be published by CMP Books in August 2001), from which this article was adapted. He enjoys receiving comments and questions from readers and can be reached by e-mail at .
2. Ledin, Jim A., “Hardware in the Loop Simulation,” Embedded Systems Programming, February 1999, p. 42.
3. Parker, T.S. and L.O. Chua. Practical Numerical Algorithms for Chaotic Systems. New York, NY: Springer-Verlag, 1989.
- Figure 1: Forward Euler integration
- Figure 2: Backward Euler integration
- Figure 3: Forward Euler integration using various values for h
- Figure 4: Example of a moderately stiff dynamic system
- Figure 5: Example of a combined discrete-continuous system
- Figure 6: Simulink diagram of turntable system with controller
- Figure 7: Digital controller model
- Figure 8: Simulation of system step response
- Figure 9: Simulated controller DAC output voltage
Return to Table of Contents