Thesis - Beaver simulation
.pdfAppendix D. Simulation of dynamic sYstems in SIMULINK- theoretical backgrounds .
D.1 Introduction.
In this appendix some theoretical aspects of the simulation routines of SIMULINKwill be outlined. The reader is advised to consult refs.[l4] and [15] for a more extensive treatment of numerical analysis. A useful book, which demonstrates the practical use of these routines specifically for aircraft simulation is ref.[29]. This reference also treats the subject of linearization and trimming of nonlinear aircraft models 2), which will not be considered here. Since the user's manual [4] does not give details about the theory behind the integration routines in SIMULINKand the software codes of the integrator routines are not accessible to the user, the equations in this appendix were reproduced from the sources mentioned above. Some differences between the expressions in this appendix and the implementation in SIMULINKmay therefore exist.
D.2 The integrators.
SIMULINKcontains six numerical integration routines: LINSIM, EULER, RK23, RK45, ADAMS, and GEAR (see also chapter 2). This section gives some general information about numerical integration methods.
D.2.1 Some fundamental concepts.
The type of problems considered.
Systems in SIMULINKneed to be described by a set of ordinary differential equations (ODEs) or ordinary difference equations:
Since few differential equations can be solved exactly, it is necessary to approximate the solutions of these ODEs numerically. Numerical integrators have been developed specifically for solving so called initid value pmblems:
This appendix is largely copied from ref.[25]. Some parts have been rewritten and some errors have been fixed.
2, |
Although trim and linearization functions are very important for flight simulation and |
|
|
|
control system design, they will not be treated in detail in this report. See refs.[ll], or |
|
[29], for more information about these subjects. |
Equation (D-2) does contain inputsignals u(t), but the methods for solving (D-1) are similar to the methods designed for solving (D-2). Problem (D-2) will be considered for the remainder of this section.
Numerical solution of ordinary differential equations (refs.[l4] and [15]).
Equation (D-2) is a vector equation. If y has N elements, N constants of integration appear in the solution. A unique solution to the system is obtained when the values of the elements of y are specified a t a certain initial time to (in this report, the time t is used as the independent variable, and to = 0). From now on, only scalar expressions will be considered. The equations can easily be transferred to more general vector equations for sets of ODES.
The SIMULINKintegrators use so called step-by-step methods, also called difference methods or discrete variable methods, to solve the set of ODES numerically. These methods generate a sequence of discrete points to, t, , ... , possibly with variable spacing h, = tn+l- t, (the stepwith). At each point t,, the solution y(t,) is approximated by y,, which is calculated from earlier values. When k earlier values yn ,Y,.~,... ,yn-k+lare used for the computation of yntl, the method is called a 'k-step method'. The Euler method (the integrator EULER in SIMULINK),for instance, is a single-step method:
Stability and errors.
Figure D-1 shows a typical family of solutions of a first order differential equation. The actual solution depends on the initial condition y(to) =yo . In figure D-1, a wrong value of yo leads to an error, which increases in time. This is called instability of the differential equation.
Initial value //
i |
Euler solution |
I |
I |
||
I |
|
I |
I |
|
I |
I |
|
I |
Figure D-1. Family of solutions of instable ODE. The Euler method is displayed graphically.
Appendix D |
169 |
In figure D-1, it is demonstrated how the Euler approximation generally crosses from one solution to another between two time steps. In this case, this leads to a n error which increases in time. In figure D-2, this growth of the error in time does not occur, because the different solutions come nearer to eachother as time proceeds. The ODE from figure D-2 is therefore stable. Nonlinear differential equations may be unstable in some regions and stable in others. For systems of (nonlinear) ODES, the situation is even more complex. One should therefore always be aware of possible instabilities of the system considered.
Figure D-2. Family of solutions of stable ODE.
Total e m
9-3.Discretization error, roundoff error, and total error as a function of h,.
There are two types of errors:
1- discretization errors (also called truncation errors),
2 - roundoff errors.
Discretization errors are a property of the numerical integration method. Roundoff errors occur due to the finite number of digits, used in the calculations, hence, they are a property of the computer and the program that is used. In general, the total error decreases a s h decreases, until a point where the round-off error becomes dominant. This is illustrated in figure D-3. Due to the errors, even for stable ODES, the numerical solution may become instable. See for instance figure D-4 and section D.2.3.
Figure D-4. Eder's method for stiff ODE.
Order of the numerical integration method.
The order of a numerical method is defined in terms of the local discretization error 6, obtained when the method is applied to problems with smooth solutions. A method is of order p if a number C exists so that:
C may depend on the derivatives of the function which defines the differen- tial equation and on the length of the interval over which the solution is sought, but it should be independent of the step number n and the step size h,. Unequality (D-4) may be written in 'big Oh notation' as:
Appendix D |
171 |
D.2.2 A review of several numerical integration methods.
It is possible to distinguish between four general categories of step-by-step methods [14], which will briefly be discussed here. All results can easily be converted to vector notations for sets of ODEs.
a. Taylor series methods.
A smooth solution y(t) can be approximated by the Taylor series expansion:
If it is possible to calculate higher-order derivatives of y(t), a method of order p can be obtained by using:
The first neglected term provides an estimate of the local discretization error. The Euler method is an example of a Taylor series method, where all derivatives of order two and higher have been neglected. Hence, the local discretization error of the Euler method is 0 (h:). Higher-order Taylor series methods involve increasingly complex computations of the higher order derivatives. For this reason, the general applicability of these methods is quite limited [14].
It is not recommended to use the EULER integrator, because it needs much smaller step sizes than the other routines in order to obtain the same accuracy. It is included in SIMULINKmainly for historical reasons [4]: it forms the base of every textbook on ODEs, because it can easily be illustrated graphically, see figure D-1. The method contained in SIMULINKuses variable step sizes, but it does not take steps back in time if the userdefined error tolerance is exceeded. The method can be converted to a fixed step size method by setting the minimum step size equal to the maximum step size.
b. Runge-Kutta methods.
Runge-Kutta methods approximate Taylor series methods, without evaluating derivatives beyond the first. However, several evaluations of the function f a r e needed to replace the higher-order derivatives. Modern RungeKutta algorithms often include techniques for estimating discretization error for step size control. The Runge-Kutta methods are self-starting, which means that only one value yn is needed to calculate Y,+~.
SIMULINKcontains two Runge-Kutta integrators, called RK23 and RK45.
RK45 is a fifth order Runge-Kutta method, which uses a fourth order method for step size control. It takes six unevenly spaced points be-
tween two output points. Four of these function values, combined with one set of coefficients, are used to produce the 4" order method. All six function values, combined with another set of coefficients, produces the 5" order method. Comparison of the two values yields the error estimate for the step size control.
RK23 is a third order Runge-Kutta method, which uses a second order method for step size control. It needs three function evaluations to generate one output point. Two steps are needed for calculating the 2ndorder method and all three function values are used for calculating the 3rd order method. The step size is again controlled by comparing the estimated error with the user-specified error tolerance.
Derivation of a second order Runge-Kutta method.
A second order Runge-Kutta method for scalar ODEs will now be investigated in detail. The results can easily be transferred to vector notation for sets of ODEs, and higher order methods can be derived in a similar way. The second order method uses two function evaluations per step. The first evaluation is:
Then a fractional step based on k l is taken:
where a and f3 are two as yet unknown coefficients. The two function values are combined to make a complete step, involving two more unknown coefficients:
To determine the coefficients, k1 and k2 are expanded in a Taylor series about (t, ,yn), giving:
Appendix D |
173 |
The expansion of yn+lcan be compared with the expression for the true local solution zn(t):
zn (tn+l ) = Z, (t,) + hn2,' (t,) + -h,2 znB'(t,) + ...
2
Matching coefficients of powers hn leads to three equations involving the unknown coefficients:
It is now possible to choose one coefficient as a parameter to produce a oneparameter family of Runge-Kutta methods. If a is chosen as a parameter this leads to:
Two obvious choises of a are % (which yields a method closely related to the rectangle rule) and 1 (yielding a method related to the trapezoid rule). The method is second order, because no choice of a can eliminate the h3 terms, neglected above (sop = 3).
The integrator RK45.
In a similar way, higher order Runge-Kutta methods can be derived. The fifth order method is described by the following formulas:
There are 27 coefficients: 6 a's, 15 p's, and 6 y's. The p's form a lower triangular array, so that each ki is obtained from previous k's. By expanding all k's in Taylor series, substituting the expansions into the formula for
Y n + l , and comparing the result with the Taylor series for ~,(t,+~),the coefficients can be determined.
For some combinations of coefficients, the difference between two expansions starts with a term including hn6,but not hn5. The resulting method is consequently fifth order @ = 6). The combination used in the Runge-Kutta Fehlberg routine are given in table D-1, see refs.[l4] and [15]. For p = 1, 2, 3, and 4, a pth order method with p function evaluations can be found. For p = 5 or p = 6, p function evaluations will produce a method of order @ - 1) only. The 'additional' function evaluation to obtain a 5th order method is not wasted, however, because a set of coefficients y i* exists, with which four of
the six k's can be combined to give a second value |
which is accurate to |
fourth order. Hence: |
|
These coefficients are also given in table D-1. For step size control the error
6
estimate C (yi - y f ) k i is used. The MATLABroutine ODE45, for example,
i=1
uses the Runge-Kutta Fehlberg algorithm, in which the estimated error 6 is compared with the acceptable error t. The estimated error 6 is computed with the vector equivalents of equations (D-15) and (D-16). The acceptable error is calculated with:
where to1 is the desired accuracy. If y is a vector, equation (D-17) becomes:
where 1y 1 is a vector in which all elements equal the absolute value of the corresponding elements in y. Knowing 6 and t, ODE45 updates the step size as follows:
ODE45 only updates the step size if 6, L 0. It is to be expected that the SIMULINKintegrator RK45 is based upon a similar algorithm.
Appendix D |
175 |
Bii
Table D-1. Fehlberg coefficients.
The integrator RK23.
The SIMULINKintegrator RK23 can be derived in a similar way. The software source code is not accessible for the user, but the algorithm is probably the same as the one that is used in the MATLABfunction ODE23. This routine uses the following expressions to obtain an approximation of y(t,):
ODE23 updates y if the estimated error 6 is smaller than or equal to the acceptable error T. These are estimated with the following expressions:
or if y and k are vectors:
The estimate of 6 involves another Runge-Kutta method than used for the update Y,+~,just as was the case with RK45.The step size is chosen as large as possible, within the user-specified limits. In ODE23, the step size is updated with the following formula:
SIMULINKuses a variable step size, but it is possible to fix the step size by setting the maximum step size equal to the minimum step size. The RungeKutta methods usually outperform the other methods when the system is highly nonlinear or discontinuous. They perform well for mixed continuous and discrete time systems. For systems with a very wide spread of time constants ('stiff systems'), LINSIM or GEAR should be used (refs.[4], [29]).
c. Multistep methods.
The methods considered so far calculate the value Y,+~ by means of a func-
tion which depends only on t y, ,and the step size h, . Multistep methods |
|
'. |
,y, |
use information a t previous points to obtain more accuracy, namely y,l |
|
,yw3, ... , and fnTl ,fw2 ,fw3 , ... , where fi = f(yi , ti). Multistep methods can be very effective. They usually require less function evaluations than one-step methods to obtain a high accuracy. Furthermore, an estimate of the discretization error can often be trivially obtained [14].All linear multistep methods can be considered as special cases of the formula:
where k is a fixed integer and either ak or Pk is not zero. The formula defines the general k-step method. It is linear, because every fi enters linearly; f doesn't necessarily have to be a linear function in its arguments.
After the method is 'started', each step requires the calculation of y,+, from the known values: y,, y,, , ... , y w k + l f,, f w l , If PO = 0, this
method is explicit and the calculation is straightforward. If Po # 0, the method is implicit, because fn+lis needed to solve for yntl.
Usually, two multistep methods are used together in computing each step of a solution: an explicit method, called predictor, is followed by one or more applications of an implicit method, which is called a corrector. These