1
Introduction

1.1 Scope of the Book

This book is for advanced undergraduate students, post-graduate students, or engineers to acquire programming skills for dynamic system modelling and analysis using control theory. The readers are assumed to have a basic understanding of computer programming, ordinary differential equations (ODE), vector calculus, and probability.

Most engineering curricula at the undergraduate level include only an elementary-level programming course in the early of the undergraduate years. Only a handful of self-motivated engineering students acquire advanced level programming skills mainly from self-study through tedious time-consuming practices and trivial mistakes. As modern engineering systems such as aircraft, satellite, automobile, or autonomous robots are implemented through inseparable tight integration of hardware systems and software algorithms, the demand for engineers having fluent skills in dynamic system modelling and algorithm design is increasing. In addition, the emergence of interdisciplinary areas merging the experimental domain with mathematical and computational approaches such as systems biology, synthetic biology, or computational neuroscience further increases the necessity of the engineers who understand dynamics and are capable of computational implementations of dynamic models.

This book aims to fill the gap in learning practical dynamic modelling, simulation, and analysis skills in aerospace engineering, robotics, and biology. Learning programming in the engineering or biology domain requires not only domain knowledge but also a robust conceptual understanding of algorithm design and implementation. It is not, of course, the skills to learn in 14 days or less as many online courses claim. To be confident in dynamic system modelling and analysis takes more than several years of practice and dedication. This book provides the starting point of the long journey for the readers to equip and prepare better for real engineering and scientific problems.

1.2 Motivation Examples

1.2.1 Free-Falling Object

Newton's second law of motion is given by

where is the -th external force in Newtons (N) acting on the object characterized by the mass, , in kg, is the time derivative, is the time in seconds, is the velocity in m/s, and is the momentum of the object. Newton's second law states that the sum of all external forces is equal to the momentum change per unit of time.

Consider a free-falling object shown in Figure 1.1. There exists only one external force, i.e. the gravitational force acting downwards in the figure. Hence, the left-hand side of (1.1) is simply given by , where is the gravitational force. Introduce the additional assumption that the object is within the reasonable range from the sea level. With the assumption, the gravitational force, , is known to be proportional to the mass, and the proportional constant is the gravitational acceleration constant, , which is equal to 9.81 m/ in the sea level. Therefore, . Replace the left-hand side of (1.1), i.e. , by provides

where the downward direction is set to the positive direction, which is the opposite of the usual convention. It highlights that establishing a consistent coordinate system at the beginning of modelling is vital in dynamic system simulation.

Figure 1.1 Free-falling object.

From the kinematic relationship between the velocity, , and the displacement, , we have

where the origin of is at the initial position of the object, , and the positive direction of is downwards in the figure. The right-hand side of (1.2) becomes

Finally, the leftmost and the rightmost terms are equal to each other as follows:

and it is expanded as follows:

Using the short notations, , , and , and after rearrangements, the governing equation is given by

For purely educational purposes, assume that the mass change rate is given by

We can identify now that there are three independent time-varying states, which are the position, , the velocity, , and the mass, . All the other time-varying states, for example, and , can be expressed using the independent state variables. Define the state variables as follows:

Obtain the time derivative of each state expressed in the state variable as follows:

and this is called the state-space form.

Let the initial conditions be equal to m, m/s, and kg. Equation (1.5) can be written in a compact form using the matrix-vector notations. Define the state vector, , as follows:

and the corresponding state-space form is written as

The second-order differential equation, (1.3), and the first-order differential equation, (1.4), are combined into the first-order three-dimensional vector differential equation, (1.6). Any higher order differential equations can be transformed into the first-order multi-dimensional vector differential equation, . Numerical integration methods such as Runge-Kutta integration (Press et al., 2007) solves the first-order ODE. They can solve any high-order differential equations by transforming them into the corresponding first-order multi-dimensional differential equation.

1.2.1.1 First Program in Matlab

We are ready to solve (1.6) with the initial condition equal to , where the superscript is the transpose of the vector. We solve the differential equation from to seconds using Matlab. Matlab includes many numerical functions and libraries to be used for dynamic simulation and analysis. A numerical integrator is one of the functions already implemented in Matlab. Hence, the only task we have to do for solving the differential equation is to learn how to use the existing functions and libraries in Matlab. The complete programme to solve the free-falling object problem is given in Program 1.1. Producing Figure 1.2 is left as an exercise in Exercise 1.1.

Figure 1.2 Free-falling object position, velocity, and mass time histories.

Now, we study the first program line by line. The m-script starts with the command 'clear'. The clear command removes all variables in the workspace. In the workspace, there would be some variables defined and used in previous activities. They may have the same names but different meanings and values in the current calculation. For example, the gravitational acceleration 'grv_const' in the third line is undefined in the current program and uses a variable of the same name used to analyse objects falling on the moon. A falling object program in the Moon was executed earlier, and 'grv_const' is still in the workspace. Without the clear command, the incorrect constant is used in the program producing wrong results. Hence, it is recommended to clear the workspace before starting new calculations. We must be careful, however, that the clear command erases all variables in the workspace. Before the clear command, we check if all values, which might be generated from a long computer simulation, were saved.

From line 3 to line 12, several constants are defined. Based on the equations we have seen earlier, it is tempting to write a code as follows:

These seem to look compact and closer to the equations we derived. It is a bad habit to write a program in this way. The list of problems in the above programming style is as follows:

• It defines a variable with a single character, 'g', 'x', 'v', etc. Using a single character variable might cause confusion on the meaning of the variable and lead to using them in wrong places with incorrect interpretations.
• Numerical numbers are written without units. There is no indication of units of the numerical values, e.g. 9.81, is it m/ or ft/?
• It uses magic numbers. What do the numbers, 0 and 5, mean in defining 't'?

Program 1.1 uses a better style. The initial position is defined using the variable name, 'init_pos', whose value is 0.0 and the unit is in metres. Appropriately named variables reduce mistakes and confusion in the program. Program 1.1 indicates the corresponding unit for each numerical value, e.g. the 'init mass' value 5.0 is in kg. We understand the meaning of each variable by its name. The texts after '%' are the comments, where we could add various information such as the unit of each numerical value.

In line 13, the built-in Runge-Kutta integrator, ode45(), is used to integrate the differential equation provided by the function, 'free_falling_obj', at the end of the m-script. Frequently, each function is saved as a separate m-script. It could also be included in the m-script for the cases that the functions might be used in the specific...