Newsletters - MATLAB Digest

Nonlinear Systems Modeling in Simulink

by Ender St. John-Olcayto

Introduction
This article shows you how to simulate nonlinear, continuous-time systems in Simulink. After describing the principles of dynamic systems, an example is given. This example is a simple continuous-time, nonlinear system and it is shown how Simulink can be used in different ways to model it. A background to the theory of system linearization is given and is applied to the example model.

Several Simulink systems have been created for use in conjunction with this article. These systems have been placed within the Simulink library sublib_digest.mdl. In addition, other MATLAB support files have been provided. These should be placed on your MATLAB's search path. All of these files can be found in this winzip archive.

The article can be read as an integrated whole or each section can be read separately. The different sections are:

A Brief Introduction to Dynamic Systems:
A background to dynamic systems is given and the state-space approach to modeling is introduced.

A Simple, Continuous-Time Dynamic System:
A simple dynamic system is described and is represented in a form suitable for modeling within Simulink

Three Ways to Model A Dynamic System:
Three main ways to model a dynamic system in Simulink are illustrated using a simple, continuous-time system. These methods are:
   1. Modeling using the standard Simulink library,
   2. Modeling using M-Code S-Functions,
   3. Modeling using C-MEX S-Functions.
A brief introduction to the use of S-Functions is also given.

Linearizing Simulink Models:
The main idea behind system linearization and linear, time-invariant systems is given and is illustrated by a simple Simulink example.

A Brief Introduction to Dynamic Systems
A dynamic system is one that evolves through time. Most systems in nature are dynamic, even if some of those appear not to change. For example, to human senses, a mountain may appear not to do anything at all! However, it is subject to wind, water and ice erosion and plate tectonics. So, its dimensions will change through time, although the time scale is measured in eons!

Another example is a thrown football. Consider the position p(t) of the football's center of mass after it has been thrown (see Figure 1). The position is dependent on several things such as its initial position and initial velocity (p₀(t) and respectively), the current velocity , the drag force acted upon it by the air (which in turn is dependent on the linear velocity) and the rotational velocity . The angle that the velocity vector makes with the relative airflow is . Additionally, the physical properties of the football, such as its mass, cross sectional area, drag coefficient and mass moment of inertia) are important. Note that the drag force F_d(t) is dependent on the properties of the atmosphere (the density) and the dynamics of the ball. That is, the system of interest is the combination of the football and the atmosphere! The system context is important: other systems involving the football include the dynamics of the football deflating though time! These dynamics may have some bearing on the flight of the football (since the football's mass will change as it deflates). However, they are likely to be insignificant for the duration of the ball's flight (unless the quarterback has a really tremendous throw!).

The only variable of interest to the spectators is p(t). Therefore, let p(t) be the output of the system. p(t) needs to be calculated from the other variables. This is often achieved by defining a state variable x(t) that summarizes the internal dynamics of a system. For this example, let x(t) be a vector-valued function defined as follows:

where the T subscript indicates the transpose operator. Furthermore, let u(t) be the force imparted by the quarterback. It is common to calculate x(t) by integrating its first-order time-derivative . is obtained by from the physics of the problem and is expressed as:

where f is a nonlinear vector function of the state x(t) and input u(t) (that is, for each element in there is a different function of x(t) and u(t)). Sometimes a change of variables is used and the state variable no longer represents physical quantites. This is sometimes done to ensure better scaling of the elements of x. Ultimately however, x is derived from the physics of the problem. The above equation is known as the state equation. The output is obtained from the following equation:

where I indicates the identity matrix and 0 is a zero matrix with dimensions indicated by the subscripts. The identity matrix picks out from x(t) the variables of interest (that is, those contained within p(t)). However, in general, the output is a nonlinear function of the state and the input:

where g, like f is a nonlinear vector function of the state x(t) and input u(t). The above is known as the output equation.

The output y is found by integrating the state equation with respect to time (the explicit dependence on time has been omitted for clarity):

and then applying the output equation. The above yields a non-unique solution due to the constant of ontegration. Thes constant is the initial state x₀ and is given by:

The actual dynamics of the football are not given here (they are quite complex!). However, the above ideas are used in the next section to model a simple nonlinear, dynamical system with a single state.

A Simple, Continuous-Time Dynamic System
Figure 1 shows a simple model of a submarine. Let it have mass M and cross-sectional area A. The operating medium is water and has a density at the operating depth and temperature. The dynamics in the x-direction is given by summing the applied forces and applying Newton's 2nd Law. The submarine is moving with a velocity V and acceleration A. Furthermore, the initial velocity of the submarine is V₀.

The drag force D is a function of the vehicle geometry and water properties. D is given by:

The applied force F is provided by the propeller. Summing forces and equating the result to the product of the vehicle mass and its acceleration gives:

A change of variables is made in order to express the above in the more general state-space representation. Let the state variable x be the velocity V (thus x₀ represents the initial velocity V₀). Using this state variable, the acceleration A is given by . Let the velocity be the measured variable y. Also, let F be written as u, the system input. The above equation can be expressed in state-space form as:

These equations can be realized in many different forms. The next section shows how these various methods can implemented in Simulink.

Three Ways to Model A Dynamic System
Introduction
Simulink is very flexible with respect to the number of different ways that a dynamic system can be represented. Dynamic systems can be represented in state-space, transfer-function and pole-zero forms. In addition, these representations can be continuous-time and discrete-time. Typically, nonlinear systems are represented using state-space equations. These equations can be realized using three different methods. This section shows how a simple continuous-time dynamic system can be represented using these methods.

Consider the simple nonlinear dynamics of the submarine model described in section 3

It can be seen that this is a first-order equation and in order for Simulink to calculate the single output y, the single state x (the vehicle velocity) must be calculated. the velocity x is calculated via:

The above equation indicates the number of integrator blocks we need. Since there is only one integration symbol, only one integrator block is required. Note the equation is implicit with respect to x. Simulink solves the above equation numerically using an integration algorithm called a solver.

Modeling Using Simulink Library Blocks
The MATLAB/Simulink environment is extremely flexible and so it comes as no surprise that there are a number of different ways of to represent systems using the block library.

A Brief Background to S-Functions
S-Functions are special MATLAB functions that permit the system dynamics to be expressed as equations. Even though Simulink blocks allow the developer to visualize the system dynamics in pictorial form, it is sometimes advantageous to write the equations of motion directly. These equations are implemented in MATLAB readable file. Within this file exists a set of functions which are known collectively as the "S-Function".

Simulink interrogates the S-Function at different phases of the integration algorithm to in order to calculate properly the outputs of the system described within it. The user must provide the means for Simulink to obtain the information that is required by it during each of these phases. Examples of how this is done will be given in the next two sections.

Since S-Functions are really MATLAB functions (which must be created with the aid of some text editor), some interface to the Simulink diagram must be provided. This interface is the S-Function block and can be found in the Nonlinear library. This block is shown in Figure 5, together with its dialog box.

The first field in the dialog box is the name of the user-written S-Function. For example, if the S-Function is implemented within the file sfun_myfunc.m then sfun_myfunc would be the S-Function name. The second field is a comma-separated list of user-defined parameters. For the submarine example, this list are the model parameters (M, A, C_d, and V₀).

When the S-Function is implemented in MATLAB's native M-Code language, these parameters are passed to the S-Function (in addition to some other right-hand side parameters that Simulink passes) just like the familiar right-hand side arguments for normal MATLAB functions. When the S-Function is implemented with a C language file (which is called a "C-MEX S-Function"), then the parameters are handled slightly differently. The next two sections will elaborate on how the parameters, and other information is passed between the S-Function and Simulink.

Modeling Using an M-Code S-Function
Using M-Code to implement is a good way of producing a prototype of the nonlinear system model. This is because the user does not need to be concerned with the complexities of the C programming language since M-Code is much easier to debug (especially with the integrated editor/debugger which works even when the M-Code function is executing within a Simulink model). Furthermore, with the advent of sub-functions in MATLAB v5.0, M-Code S-Functions can be made to look like their C-MEX counterparts.

The S-Function sfun_sub.m (in file sfun_sub.m) implements the submarine dynamics. Open it using the editor and inspect the different sections. The first line of the function shows how the user-defined parameters are passed to the function. They appear as additional right-hand side arguments following the t, x, u, flag arguments. The latter are passed to sfun_su by Simulink during the different phases of the integration algorithm. These arguments correspond tob the current time, state, input and a simulation status flag. It should be noted that the current time does not necessarily correspond to the simulation clock time (which can be viewed using the clock block in the Sources library). This is because each integration time-step (called the "major time-step" is performed using a series of smaller steps, each of which is called a "minor time-step").

Simulink interrogates the S-Function during the different integration phases by setting the value of the flag argument. This argument determines whether Simulink requires information on the initial conditions, the continuous-time state derivative, the next state (for discrete-time systems), the output and so on. The user must test for the value of flag and supply the appropriate algorithm accordingly. For example, when flag=1, then Simulink will calculate using the user-supplied state equation. Similarly, when flag=3, then Simulink will calculate using the system's output equation. The additional arguments supplied by Simulink (t, x and u) can be used within each of these user-supplied algorithms.

Simulink uses the left-hand argument sys in order to return to the values calculated by each of the algorithms. This means that sys takes on different meanings according to the value of flag.

The x0 argument is the initial state of the system represented within the S-Function. The initial state is considered to be too important to be placed within sys! This should be specified in the initialization section: the portion of M-Code that is executed when flag=0.

str should be set to the empty matrix within the initialization section. It exists only to comply with the S-Function application programmer interface for block diagrams. Apart from the assignment str=[], this argument can be ignored!

The last left-hand side argument ts, is a matrix of sample time information and is specified by the user within the initialization section. Note that sample time information must be supplied even if the system being modeled is continuous-time. Suppose the system being modeled has three states. Let the first state x_c be continuous-time and let the other two be discrete-time (x_d1 and x_d2). Furthermore, x_d1 and x_d2 have sample times T₁ and T₂ respectively. However, the first sample of x_d1 does not occur at t=0. Rather it is offset by an amount . Similarly, the first sample of x_d2 is offset by an amount . For this case, ts is specified to be:

The next section will discuss the main points regarding the writing of S-Functions.

Modeling Using a C-MEX S-Function
Why C-MEX? C because the S-Function is written in the C programming language and MEX because a mex-file is generated when the C-file is compiled. A mex-file is the generic name for a dynamic link library that implements a MATLAB function. The generated file has the same name as the C source-code file. However, they have different filename extensions according to the platform on which they are implemented. For PC/Windows platforms the extension is dll whereas on Solaris platforms the extension is mexsol.

S-Functions come in two main flavors:
• Level 1 S-Functions
• Level 2 S-Functions

Only the C language is capable of implementing level 2 S-Functions. S-Function mex-files can be implemented in FORTRAN but more flexibility can be obtained if they are written in C. However, it is always possible to call FORTRAN routines from C!

Before discussing how to create the mex-file, the file sfun_sub.c will be examined and compared to sfun_sub.m. Things to note are:
• The filenames are the same except for the filename extension:

This permits the compiled S-Function to be called in exactly the same way as the corresponding M-Code S-Function. No changes are required to Simulink S-Function blocks: MATLAB will use the mex-file version of a function in preference to an M-Code version.
• There is some commonality between the names of the functions in sfun_sub.c and the sub-functions in sfun_sub.m:

This has been done by design! This is so that S-Functions can be prototyped quickly using MATLAB's native language.
• The name of the C-MEX S-Function is specified differently from the M-Code version:

The name of a C-MEX S-Function is specified using a C pre-processor definition such as

#define S_FUNCTION_NAME sfun_sub

This means that a C-MEX S-Function can be defined within a file with a different filename. • The C-MEX S-Function is of type level-2:

This is specified using the C pre-processor definition:

#define S_FUNCTION_NAME sfun_sub

It is recommended that all new C-MEX S-Functions be level-2. • Most of the internal functions accept only one argument:

This argument is called a pointer to a SimStruct. A SimStruct is a C language data-type that contains all the information that Simulink needs to carry out the simulation of the system implemented within the S-Function. • There are lots of C pre-processor macros that act on the SimStruct:

These macros let the user get and set information regarding the system. For example, some of these macros are used to obtain the state derivatives and outputs. Others are used to set the numbers of inputs and outputs. The SimStruct is never manipulated directly. This is because The MathWorks Inc. may need to change the underlying architecture of the SimStruct as new features are developed. In order to prevent the user from having to modify their code when these changes are made, explicit reference to the SimStruct is avoided via the use of the macros. A list of the SimStruct access macros can be found in the Using Simulink guide. • There does not appear to be any mechanism to provide left-hand side and additional left-hand side arguments:

Left-hand and right-hand side arguments do not need to be supplied since the system information is conveyed using the S-Function's SimStruct. Additional, user-defined parameters are obtained from the S-Function block parameter list via the macro ssGetSFcnParam. However, left-hand and right-hand side arguments may need to be supplied if the S-Function is executed from the MATLAB command-line. The mechanism to support this is contained within the file simulink.c. This file is included after all of the user-written code. Look for the following line in sfun_sub.c (it can be found near the end of the file):

#include "simulink.c"

The mex-file can be created by issuing the mex sfun.c command within MATLAB. Additional support files may be compiled and linked by providing them as additional arguments to the S-Function.

Masking the Model
The final task of creating a custom block using S-Functions is to create a dialog box so that the user can supply different parameters to it (for example, it may be required to simulate a submarine with a different mass).

Masking an S-Function requires three steps: • Replacing hard-coded parameters in the S-Function block with variables • Encapsulating the S-Function block within a subsystem • Providing the dialog box using the Mask Editor (found under the Edit drop-down menu of the model)

Figure 6 shows the actions taken to create the mask for the simple submarine model.

Linearizing Simulink Models
Introduction
What does linearization mean and why are nonlinear systems linearized? The first question is best answered by explaining what a nonlinear system is. Most real systems are nonlinear. The speed of an automobile does not double if the gas pedal is pressed twice as hard, or the Earth's average surface temperature will not increase at the same rate at which greenhouse gases are injected into the atmosphere. Rather, the output behavior is a complex nonlinear function of the inputs into the system (the automobile and the atmosphere in the above examples).

The outputs may also depend on the time-history of some internal variables (called states). In this case the system is a nonlinear dynamic system (the submarine model is an example of a nonlinear dynamic system with one state).

Some systems are yet more complex. Suppose that a system S is excited (or "forced") by an input u(t) at time T₁ and that the resulting output response is y₁(t). Now suppose that S is forced by the same input but at some later time T₂ (this input u(t+T₂) is the same shape as u(t) but shifted by T₂ time units). The output of S in this case is y₂(t). If y₂(t)₁(t+T₂), then S is time-varying. A familiar time-varying system is the automobile. As time progresses, the mass of the vehicle decreases due to the consumption of fuel during the journey. This means that as time progresses, the vehicle will accelerate more for a given gas pedal position. Time-varying systems will not be discussed in this article but will be presented in a future article.

Sometimes it is desired to maintain the nonlinear system about some "operating point" (sometimes called a "trim point"). A nonlinear system may have several equilibrium points. Consider another system S which has three equilibrium points E₁,E₂ and E₃. Figure 7 shows some state trajectories of S that start from the initial conditions x₀₁, x₀₂ and x₀₃ and from E₃. When the system starts from x₀₁, then the system comes to rest at equilibrium point E₁ (that is, the state derivative ). When the system starts from x₀₂ or x₀₃, the system comes to rest at E₂. Note though, that if the system starts from E₃ and is given an input u₁ then it will come to rest at E₁. Similarly, if the system starts from E₃ and is given an input u₂ then it will come to rest at E₃. Equilibrium points that attracts trajectories are called "stable" equilibrium points. Similarly, equilibrium points from which trajectories emerge are called stable equilibrium points. Note that in order for a trajectory to emerge from E₃, the system had to be forced by inputs u₁ and u₂. If there had been no input, then S would have remained at the unstable equilibrium point E₃.

Linearization means approximating the motion about each equilibrium as being linear. Consider the nonlinear, time-invariant system S at one of its equilibrium points E. The state and input at E are x₀ and u₀ respectively. The state and output equations at E are given by:

Let S be perturbed by a small input where (that is, the total input is ). This causes the state derivative and output to be perturbed resulting in:

Approximating the above using a Taylor series about x₀ and u₀ gives:

where HOT means higher order terms (those in the series beyond the first order terms). It has already been established that at E, . Considering only the perturbation outputs around the equilibrium output y₀ and by recognizing that the terms involving partial derivatives are matrices, yields the linear equations:

where A, B, C and D are given by:

It is more usual to drop the symbols since it is taken as known that the above equations represent the perturbation dynamics of S about the equilibrium E. Therefore the state and output equations become:

The above is a Linear Time-Invariant system. Note that if S is a time-varying nonlinear system ( and ), then the partial derivatives would be carried out with respect to t in addition to x and u. In this case, the matrices would be time-varying: A=A(t), B=B(t), C=C(t), D=D(t) and the linearization would yield a Linear Time-Varying system. The next section discusses how the above procedure is automated within Simulink.

Linearizing a Continuous-Time System
If a nonlinear system has been implemented in Simulink, then it can be linearized by first finding the equilibrium points and then linearizing it about those equilibrium points using the MATLAB functions trim and linmod respectively. The masked subsystem sublin_digest can be used to linearize the submarine model. Run the script trimsub_at_1 to linearize the submarine at 1 m/s. This script opens the model sublin_digest (see Figure 8) and sets the submarine parameters using the function mskvar (not a standard toolbox function!). It then trims the system at the equilibrium point corresponding to 1 m/s using trim and linearizes it using linmod. In addition, the resulting state space model is converted to two LTI objects (the Control System Toolbox is required in order to do this!).