To read original PDF of the print article, click here.
Modeling Dynamic Systems
To develop a simulation of a complex dynamic system, you must first develop mathematical models of major system components, as well as of any significant interactions between the system and its operational environment. Here's an introduction to the development of mathematical models of dynamic systems.
A mathematical model is an algorithm or set of equations that is combined with a set of data values to represent the significant behavior of a system, process, or phenomenon.
The development of a mathematical model for a given real-world system can be a difficult task. In cases where the system's dynamics are not well understood, a series of experiments must be performed to collect data that can then be processed using various techniques to yield a model of system behavior. This article introduces some of the methods used in the development of mathematical models of real world systems and phenomena.
Continuous-time dynamic systems
Some examples of continuous-time dynamic systems are the translational and rotational motion of an aircraft, the orbital motion of a satellite, the response of a robotic arm to the motion of its actuators, and the response of an op-amp bandpass filter to an input signal. The differential equations describing the behavior of these systems are called dynamic equations.
High-fidelity dynamic equations representing complex, real-world systems tend to be nonlinear and time-varying, which makes their analytical solution difficult or impossible. When simulating these systems, numerical integration algorithms are used to estimate the solution of the dynamic equations. Because numerical integration is used, a mathematical model of a dynamic system only requires the dynamic equations that describe the system's behavior. The actual solution of these equations occurs while the simulation is running.
Most numerical integration algorithms operate only on first-order differential equations. This means that if a dynamic equation contains second-order (or higher) derivatives, it must be transformed into an equivalent set of first-order differential equations. This is a simple procedure, which is described in the sidebar ≥Modeling High-Order Differential Equations.≤
Equation 1 below shows a second-order (nonlinear) differential equation. The result of solving Equation 1 for the highest derivative is shown in Equation 2. Equation 3 is a set of two first-order differential equations that are equivalent to Equation 1. In Equation 3, x1 is equal to the solution function x, and x2 is equal to the derivative x'. The variables x1 and x2 are referred to as state variables:
In general, a differential equation such as Equation 1 will have many solution functions. Additional information must be provided to specify the solution of interest. In dynamic system simulation, this additional information is given in the form of initial conditions on the state variables. These initial conditions are the values of the state variables at the beginning of the integration interval.
Equation 4 shows an example set of initial conditions at time zero that, combined with Equation 3, uniquely specify the solution to the dynamic equation of Equation 1:
The process of model development begins by specifying the requirements the model must meet. Some of the issues that must be addressed are:
These issues should be addressed as part of planning for the simulation effort. The questions above can be applied initially at the highest level of the entire system being simulated and then again and again as the system is broken down into subsystems and individual components to be modeled. The same questions should also be used in the development of additional models that are required for a complete simulation, such as the gravitational field and solar pressure models in the communication satellite previous example.
Level of model complexity
One useful approach to dealing with these issues is to start with a relatively simple model containing a limited set of effects and a coarse level of model detail. Then, as experience is gained with the simulation, more effects and model details can be added as needed. Often, in the early stages of simulation development, it is not clear which effects and model details are truly significant. If a large number of effects and model details are included in the initial design of the model, it may turn out that much effort has been wasted modeling things that turn out to be insignificant in determining system performance.
If the software interfaces to each model are well defined, it should be possible to replace individual models with higher fidelity representations without the need to make significant changes to the rest of the simulation. It may also be useful to maintain multiple levels of fidelity simultaneously for particular models in the simulation. This will allow the simulation user to select the desired level of fidelity as part of the simulation input data set. Having multiple model fidelity levels available in a simulation can allow the user to perform detailed modeling of particular effects or system details when needed, and potentially reduce execution time significantly on simulation runs when that level of detail is not required.
Many systems have behavior that is too complex to represent easily in terms of the laws of physics. The aerodynamics of a supersonic aircraft, for example, tend to be complex and nonlinear. In this case, the only reasonable approach may be to measure the system's behavior with a sub-scale model in a wind tunnel and then create a set of lookup tables to serve as the model. This is called an ≥empirical≤ model.
Let's look at an example of a physics-based model. Figure 1 shows a pendulum suspended from a string of length l under the influence of gravitational acceleration g. The pendulum angular deflection with respect to the vertical is q, given in radians. The mass of the pendulum bob is defined as m. The goal for this model will be to determine the period of oscillation of the pendulum as a function of the initial deflection angle q0 assuming that the initial velocity, q0, is zero.
Given the goal of determining the oscillation period, the effects that must be modeled will now be considered. We will look at the relevant physical effects and determine which to include in the model and which to ignore:
We know that a real pendulum will eventually slow down and stop due to friction in the string and air resistance. This isn't the kind of behavior we are interested in, so we will modify our goal to be the determination of the oscillation period at the time its motion is started. This assumes that the pendulum slows gradually and the oscillation period changes slowly.
We have made several simplifying assumptions that will make the model development task easier. Next, we will apply the laws of physics to the system to develop the dynamic equations.The force of interest due to gravity will act on the pendulum bob in the direction perpendicular to the string at any moment in time. This force is defined as shown in Equation 5:
Applying Newton's law F = ma to the problem leads to Equation 6, where a is the acceleration of the bob in the tangential direction:
The acceleration a is related to the angle q by the relation a = lq. This leads to the final dynamic equation given in Equation 7:
Note that this equation does not depend on the mass of the bob m; however, it does depend on the assumptions we listed previously. It is also a nonlinear differential equation because the term sinq appears in it. To completely determine a solution for this equation the initial conditions of the system must be specified as shown in Equation 8. The parameter q0 in Equation 8 is the angle from which the bob is released at time zero with an initial velocity of zero. For this example, the possible values for q0 will be assumed to lie in the range [≠p/2, p/2]:
To make Equation 7 suitable for simulation, it must be restated as a set of first-order differential equations and corresponding initial conditions as shown in Equation 9:
Using the techniques of numerical integration, these equations can be solved for various values of q0 and the oscillation period can be determined from examining the solutions.
This example has demonstrated the basic approach for developing a physics-based mathematical model of a dynamic system. Similar techniques can be used in the development of models from other disciplines such as electronics and chemistry, as long as the dynamic equations describing the system are well defined and the data values appearing in the equations can be determined with sufficient precision. In this example, the only data values required were the gravitational acceleration g, the string length l, and the initial displacement q0. Of course, other data would have to be examined to verify that the assumptions are reasonable, such as the assumption that the mass of the string is small relative to the mass of the bob.
In situations where the dynamic equations or data values for a mathematical model are not known to the desired degree of precision, alternative methods for model development are necessary. These techniques are discussed in the next section.
Table interpolation is a static modeling technique used to evaluate functions of the form shown in Equation 10:
This is a static method in the sense that it does not permit the direct implementation of dynamic equations. However, table interpolation functions can be used in the construction of dynamic equations. For example, coefficients appearing in the dynamic equations can be computed using table interpolation.
This approach is normally used when the output of the function must be determined experimentally. It is also applicable as a speed optimization technique if a lengthy computation (perhaps an iterative procedure) is required to evaluate the function. In this case, a table interpolation to estimate the result of the computation may execute many times faster than a direct computation.
The function inputs x1, x2, and so on can be any variable in the simulation, such as time, a state variable, or even a constant. The number of function inputs is arbitrary, but in practical applications it is usually five or less. As more inputs are added to the function, its memory requirements and execution time will tend to increase. The output y depends only on the values of these inputs when the function is evaluated.
A one-dimensional example of the data set for an interpolation table with eight equally spaced breakpoints is shown in Figure 2. The span of the input variable x is [0, 0.7]. If the input variable precisely matches the x location of one of the breakpoints it is a simple matter to return the corresponding y value as the result of the function evaluation. If the input value falls between the breakpoints an interpolation must be performed.
Many different techniques are available for interpolating functions between breakpoints. We will discuss the commonly used method of linear interpolation. One-dimensional linear interpolation is performed graphically by drawing straight lines between adjacent points as shown in Figure 3. The interpolated function changes continuously between breakpoints, but its derivative is discontinuous at the breakpoints. If this not acceptable, other interpolation techniques that produce smoother results (such as cubic spline interpolation1 ) can be used instead. The steps in linear interpolation for the case where x coordinates are equally spaced are shown in the sidebar ≥Linear Interpolation with Equally Spaced Breakpoints.≤
Unequally spaced breakpoints
Although the bisection method is the most general technique for locating the correct breakpoints interval, it may be possible to eliminate this search much of the time. If the input value x changes slowly between evaluations of the function, the simple step of checking if the input is contained in the same interval as it was for the previous function evaluation will often eliminate the need for bisection. If the input does not lie in the same interval as was used in the previous evaluation, bisection would then be performed.
The bisection algorithm is considerably more time consuming than the direct computation used with equally spaced breakpoints. The advantage of using unequally spaced breakpoints is that it may be possible to adequately model a particular function with a much smaller table than would be required with equally spaced breakpoints. This is because the unequally spaced breakpoints can be closely spaced in regions where the function has rapid fluctuations and more widely spaced in regions where the function is relatively smooth. When using equally spaced breakpoints, all the breakpoints must be spaced closely enough to accommodate the most rapid function fluctuations in the function even if this only occurs over just a small part of the input variable's span.
The selection of equally spaced vs. unequally spaced breakpoints depends on the characteristics of the function to be interpolated, the time available for function evaluation, and the memory available for table storage. If unequally spaced breakpoints are used, the locations of the breakpoints should be selected carefully to minimize the number of breakpoints required while simultaneously limiting the magnitude of the interpolation error. When equally spaced breakpoints are used, the point interval must be chosen to limit the maximum interpolation error to an acceptable value.
Multidimensional table interpolation
This example uses equally spaced breakpoints, but the steps are identical for unequally spaced breakpoints once the correct intervals have been determined. In Figure 4, the function inputs x and y define the point p, at which we wish to evaluate the function output z. The four points defined by the intersection of the lines x=x1, x=x2, y=y1, and y=y2 represent the surrounding breakpoints defined in Equations 14, 15, and 16. Two-dimensional interpolation is performed in two steps.
First, the interpolated function values at the points p1 and p2 in