Modeling 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
The behavior of a dynamic system evolves over time. This behavior is usually represented by differential equations when modeling continuous-time systems. A system is called continuous-time if its descriptive equations are defined for all values of time. Our focus in this article will be on the modeling of continuous-time 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:

Mathematical modeling
Mathematical models are developed using the techniques of engineering disciplines relevant to the system being modeled. Experts with thorough knowledge of the system and components of interest typically perform the model development.

The process of model development begins by specifying the requirements the model must meet. Some of the issues that must be addressed are:

• What effects should be included in the model? A complex system exhibits many different kinds of behavior (for example, the motion of motors, vibration, wear of moving parts, and so on), but not all of these behaviors may need to be modeled to produce an effective simulation. Limiting the effects modeled to only those that are necessary will make the model less complex and easier to build, test, and maintain, as well as requiring less computational resources to execute
• How detailed must the model be? In many cases, a simple model is all that is needed, but if precise determination of system behavior is required, the model may need to be very complex
• What interactions between the system and the outside environment must be modeled? For example, a model of communication satellite motion must typically operate in conjunction with a model of the earth's gravitational field, as well as models of other relevant phenomena such as solar pressure
• What techniques will be used to develop the model? A fundamental choice is whether to use physics-based equations or measured data as the basis for the model, or some combination of the two. The answer to this question is often obvious to those with expert knowledge of the system being modeled
• What data must be gathered to perform the modeling? For example, an aerodynamic model of an aircraft may require extensive wind tunnel testing
• How much time and how many people are available to develop the model? As model complexity increases, so do development and test hours
• What computing resources are available for the model? A large model may consume significant amounts of memory, disk space and CPU time. However, given the capabilities of current computers, this may not be a critical issue
• Will the model eventually be used in a hardware-in-the-loop (HIL) simulation? This may place severe constraints on the execution time allowed for the model. Alternatively, a complex model may require high performance computing hardware for use in a HIL simulation, perhaps involving the use of multiple processors
• How will the model implementation be verified and validated? Reasonable ways must be found to confirm that the model is correctly implemented and that its behavior matches the real-world system being modeled to an acceptable degree

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
The required complexity of a mathematical model is primarily determined by the answers to the first two questions listed in the previous section: the effects to be modeled and the level of modeling detail required. For most systems complex enough to make a simulation worthwhile, a large number of effects can be identified that potentially have some bearing on performance. The model developer must determine which effects are significant and which can be ignored. This is partially an economic decision because the more effects that are added to a model, the more it will cost to develop and validate, and the longer it will take to complete development.

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.

Modeling methods
We will now discuss some of the commonly used techniques for developing the equations and data sets for a mathematical model. A model is called ≥physics-based≤ if it is based on the equations of generally accepted physical laws. A spacecraft orbital model based on Newton's laws of motion is an example of a physics-based model.

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:

• Gravity must be modeled, since the pendulum would not move at all without it
• The mass of the pendulum bob must be modeled for the same reason
• If the size of the bob is small in comparison to the length of the string, the bob can be modeled as a point mass. This simplifies the model significantly
• If the mass of the string is much less than that of the bob, the string's mass can be ignored
• Friction within the string will be assumed to be a small effect over short time periods and will be ignored
• The pendulum will be assumed to move slowly so that air resistance is not a significant factor over a short time period

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.

Empirical modeling
Empirical modeling techniques use measured data from various types of experiments to develop a mathematical model of a system. In reality, all mathematical models are empirical to some degree. For example, the pendulum model in the previous section includes some empirically determined constants. However, our interest in this section is on the development of models for systems with substantial dynamic behavior that is not readily modeled by well-defined equations.

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 interpolation¹ ) 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
If the x coordinates of the breakpoints are not equally spaced, it takes more work to determine which point interval contains the function input value. A general approach for locating the correct interval is the technique of bisection. This is an algorithm for performing an efficient search of an ordered list. The steps in the bisection method are shown in the sidebar ≥Linear Interpolation with 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
We will now examine linear interpolation of multiple-input functions. To evaluate a multiple-input interpolation function, each input variable must have a set of breakpoints (either equally or unequally spaced) associated with it, as well as a table of the function values at the breakpoints. There is no reason that a particular function cannot have equally spaced points for some inputs and unequally spaced points for others. For each input, the interval containing the current input must be located using the interpolation techniques described previously. Then, using these breakpoints, a multidimensional interpolation is performed. We will look at an example using a two dimensional table with equally spaced breakpoints for both input variables. It is straightforward to extend this example to three or more dimensions. Equations 14, 15, and 16 list the x and y axis breakpoints and the function values z at those breakpoints. The table in Equation 16 is defined so that increasing values of x appear across the columns from left to right and increasing values of y lie along the rows from top to bottom:

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