INDEX

Tutorial Name of the Tutorial Page No.

No.

1 Introduction to Control Systems in MATLAB/Simulink 09

2 Introduction to Simulink 15

3 Solve the systems of linear equations 46

4 Modeling and Simulation of Mass-Spring-Damper System Using

47
Simulink

5 Interaction of a Simulink Model from a MATLAB Script 62

6 Modeling and Simulation of a DC Motor Using Simulink 68

7 Design and Implementation of PID Controller in Simulink 69

8 Extraction of a linear Model into MATLAB 88

9 Implementation of Kalman Filter in Simulink 104

10 Modeling and Simulation of Aerospace Vehicles 110

 Tutorial 1: Introduction to Control Systems in MATLAB/Simulink

Example 1.1: Generate a Pole-Zero map of dynamic systems with customizable plotting
options.

H (s) 2s2 + 5s + 1 1s2 + 4s + 1

G (s ) = = and 2
F ( s ) 1s 2 + 3s + 5 4s + 2s + 5

Flowchart:

MATLAB Code:

clc;
clear;
system1 = tf([2 5 1],[1 3 5]);
system2 = tf([1 4 1],[4 2 5]);
plotoptions = pzoptions;
plotoptions.Grid = 'on';
pzplot(system1,'r',system2,'g',plotoptions);
a = findobj(gca,'type','line');
set(a(1),'markersize',10)
set(a(2),'markersize',10)
set(a(3),'markersize',10)
set(a(4),'markersize',10)
set(a(5),'markersize',10)
ax = gca;
li = ax.Children(1).Children;
for k = 1:numel(li)
% li(k).LineStyle = "-";
li(k).LineWidth = 2;

2
end
li = ax.Children(2).Children;
for k = 1:numel(li)
% li(k).LineStyle = "-";
li(k).LineWidth = 2;
end
print('Example_3_1','-dpng','-r300');

Results:
Simulation of two dynamic systems with different color indications.

Example 1.2: Plot the root locus in MATLAB for a transfer function with a variable.

Case I:

K
G (s ) H (s ) =
s ( s + 3 ) ( s + 3s + 4.5 )
2

Case II:

K
G (s ) H (s ) =
s ( s + 2 ) ( s 2 + 8s + 20 )

3
Flowchart:

MATLAB Code:

%% Case_I

clc;
clear all;
k = 0:0.5:20000;
a = 1;
A = 1;
B = 6;
C = 13.5;
D = 13.5;
E = 0;
sys1 = tf(a,[A B C D E]);
h = rlocusplot(sys1,'r',k);
p = getoptions(h); % Get options for plot.
% p.Title.String = 'Root-locus Plot'; % Change title
setoptions(h,p); % Apply options to plot.
a = findobj(gca,'type','line');
set(a(1),'markersize',10)
ax = gca;
li = ax.Children(1).Children;
for k = 1:numel(li)
% li(k).LineStyle = "-";
li(k).LineWidth = 2;
end
4
print('Example_3_2','-dpng','-r300');

Root locus plot of a dynamics system with variables in the Transfer Function:

%% Case_II

clc;
clear all;
k = 0:0.5:20000;
a = 1;
A = 1;
B = 10;
C = 36;
D = 40;
E = 0;
sys1 = tf(a,[A B C D E]);
h = rlocusplot(sys1,'r',k);
p = getoptions(h); % Get options for plot.
% p.Title.String = 'Root-locus Plot'; % Change title
setoptions(h,p); % Apply options to plot.
a = findobj(gca,'type','line');
set(a(1),'markersize',10)
ax = gca;
li = ax.Children(1).Children;
for k = 1:numel(li)
% li(k).LineStyle = "-";
li(k).LineWidth = 2;
end

5
print('Example_3_2','-dpng','-r300');

Root locus plot of a dynamics system with variables in the Transfer Function:

Example 1.3: Create a Bode plot using a transfer function in MATLAB.

100 (10s + 1)
Case I : G ( s ) H ( s ) =
s ( s + 0.4 )( s + 10 )( s + 1)

4 ( 0.5s + 1)
Case II : G ( s ) H ( s ) =
s ( 2s + 1) ( 0.0156s 2 + 0.05s + 1)
2

MATLAB Code:

clc;
clear all;
sys1 = tf([1000 100],[1 11.4 14.4 4 0]);
sys2 = tf([2 4],[0.03125 0.1156 2.05 1 0 0]);
bode(sys1,'b',sys2,'r:');
legend('System_1','System_2')
grid;
s = allmargin(sys1);
t = allmargin(sys2);
h = findobj(gcf,'type','line');
set(h,'linewidth',2)

print('Example_3_3','-dpng','-r300');

6
7
 Tutorial 2: Introduction to Simulink
Simulink® is a graphical modeling and simulation environment for dynamic systems. You can
create block diagrams, where blocks represent parts of a system. A block can represent a
physical component, a small system, or a function. An input/output relationship fully
characterizes a block.

Open Simulink from the MATLAB command window:

In the MATLAB® Command Window, Type simulink to open the platform.

>> simulink

Fig. 2.1 MATLAB Page with Command Window

From MATLAB Home Page:

Open MATLAB > GOTO HOME Tab > New – Simulink Model

Fig. 2.2 MATLAB Home Page

8
Alternatively, you can hit the Simulink button at the top of the MATLAB window as
shown here:

Fig. 2.3 Home Page of MATLAB

When it starts, Simulink brings up a single window, entitled Simulink Start Page which can
be seen here.

Fig. 2.4 Page to create Blank Model

Once you click on Blank Model, a new window will appear as shown below.

Fig. 2.5 Simulation Page of SIMULINK

9
In Simulink, a model is a collection of blocks that, in general, represents a system. In addition
to creating a model from scratch, previously saved model files can be loaded either from
the File menu or from the MATLAB command prompt.

Simple.slx

Open this file (the file will be shared during lab hours) in Simulink by entering the following
command in the MATLAB command window. (Alternatively, you can load this file using
the Open option in the File menu in Simulink, or by hitting Ctrl-O in Simulink).

>> simple

The following model window should appear.

Fig. 2.6 Simulation Page with Simulink Model

A new model can be created by selecting New from the File menu in any Simulink window (or
by hitting Ctrl-N).
Basic Elements:
There are two major classes of items in Simulink: blocks and lines. Blocks are used to
generate, modify, combine, output, and display signals. Lines are used to transfer signals from
one block to another.

Blocks:

There are several general classes of blocks within the Simulink library:

▪ Sources: used to generate various signals

10
▪ Sinks: used to output or display signals
▪ Continuous: continuous-time system elements (transfer functions, state-space models,
PID controllers, etc.)
▪ Discrete: linear, discrete-time system elements (discrete transfer functions, discrete state-
space models, etc.)
▪ Math Operations: contains many common math operations (gain, sum, product, absolute
value, etc.)
▪ Ports & Subsystems: contains useful blocks to build a system

Blocks have zero to several input terminals and zero to several output terminals. Unused input
terminals are indicated by a small open triangle. Unused output terminals are indicated by a
small triangular point. The block shown below has an unused input terminal on the left and an
unused output terminal on the right.

Fig. 2.7 Transfer Function Block

Lines:

Lines transmit signals in the direction indicated by the arrow. Lines must always transmit
signals from the output terminal of one block to the input terminal of another block. One
exception to this is a line can tap off of another line, splitting the signal to each of two
destination blocks, as shown below.

Fig. 2.8 Simulink Model with Connecting Lines

Lines can never inject a signal into another line; lines must be combined through the use of a
block such as a summing junction.

11
A signal can be either a scalar signal or a vector signal. For Single-Input, Single-Output (SISO)
systems, scalar signals are generally used. For Multi-Input, Multi-Output (MIMO) systems,
vector signals are often used, consisting of two or more scalar signals. The lines used to
transmit scalar and vector signals are identical. The type of signal carried by a line is
determined by the blocks on either end of the line.

Fig. 2.9 Simple Model with Step, Transfer function, and Scope

The simple model consists of three blocks: Step, Transfer Function, and Scope. The Step is
a Source block from which a step input signal originates. This signal is transferred through
the line in the direction indicated by the arrow to the Transfer Function Continuous block.
The Transfer Function block modifies its input signal and outputs a new signal on a line to
the Scope. The Scope is a Sink block used to display a signal much like an oscilloscope.

Modifying Blocks:

A block can be modified by double-clicking on it. For example, if you double-click on

the Transfer Function block in the Simple model, you will see the following dialog box.

12
 Fig. 2.10 Block Parameters for Transfer Function

Simulink provides block libraries that are collections of blocks grouped by functionality.

Fig. 2.11 Library Browser to select the required blocks

Gathering the Blocks:

Follow the steps below to collect the necessary blocks:
▪ Create a new model (New from the File menu or hit Ctrl-N). You will get a blank model
window.
▪ Click on the Tools tab and then select Library Browser.
▪ Then click on the Sources listing in the Simulink library browser.

13
 Fig. 2.12 Simulink Library Browser

Commonly Used Blocks Library:

SIMULATION > Library Browser > Simulink > Commonly Used Blocks

For example, you use a Gain block from the commonly used blocks to model a system that
multiplies its input by a constant.

Fig. 2.13 Gain Block

To model the sine wave input to the system in Simulink, include a Sine Wave source.

Fig. 2.14 Sine wave Block (Source) + Gain Block

14
 Fig. 2.15 Commonly Used Blocks in Simulink

Inport blocks are ports that serve as links from outside a system into the system. Outport
blocks are output ports for a subsystem. A Subsystem block represents a subsystem of the
system that contains it.

Example 2.1: Modeling the differential equation with Simulink

d2 x dx
2
= −4 − 3x + 3u0 ( t )
dt dt
Required block(s):
• Step (input) and Scope (output)
• Gain
• Integrator
• Sum

Fig. 2.16 Required Blocks

15
Block diagram for the equation:

Fig. 2.17 Block Diagram for Example 2.1

Simulink Model:

Fig. 2.18 Model for Example 2.1

As a first step, we enclose the blocks and connecting lines that we want to include in the
subsystem within a bounding box. This is done by forming a rectangle around the blocks and
the connecting lines to select them. Keep the cursor on three dots at the right-down corner
and select Create Subsystem. Simulink replaces the selected blocks and connecting lines with
a Subsystem block as shown in Figure below.

Fig. 2.19 Blocks and Lines Selected to Create Subsystem

16
 Fig. 2.20 Model with Subsystem Block

Double-click on the Subsystem block in Figure 2.20, and we observe that Simulink displays
all blocks and interconnecting lines as shown in Figure 2.19 where the Step and Scope blocks
in Figure 2.19, have been replaced by In1 and Out1 blocks respectively.

Fig. 2.21 Model with Inport and Outport Blocks

1.1 Ground Block:

The Ground block can be used to connect blocks whose input ports are not connected to
other blocks.

Example 2: Perform the operation A + B given A = 3+2j and B = 4+3j

Fig. 2.22 Display of the sum of two complex numbers

17
If we run a simulation with blocks having unconnected input ports, Simulink issues warning
messages. We can avoid the warning messages by using Ground blocks.

Fig. 2.23 The model with Input Block Deleted

Fig. 2.24 Model with a Ground block connected to the Sum block

1.2 Terminator Block:

The Terminator block can be used to cap blocks whose output ports are not connected to
other blocks. If we run a simulation with blocks having unconnected output ports, Simulink
issues warning messages. We can avoid the warning messages by using Terminator blocks.

Fig. 2.25 Sum Block with Unconnected Output

18
 Fig. 2.26 Sum block with output connected to a Terminator block

1.3 Constant and Product Block:

Constant Block:

The Constant block is used to define a real or complex constant value. This block accepts
scalar (1x1 2−D array), vector (1−D array), or matrix (2−D array) output, depending on the
dimensionality of the Constant value parameter that we specify, and the setting of the Interpret
vector parameters as a 1−D parameter.
Product Block:

The Product block performs multiplication or division of its inputs. This block produces
outputs using either element−wise or matrix multiplication, depending on the value of the
Multiplication parameter.

Example 2.3: Perform the Multiplication ( 3 + j4 )  ( 4 + j3 )  ( 5 − j8 )

Fig. 2.27 Model with Product Block

19
1.4 Divide Block:

Example 2.4: Perform the following operation ( 3 + j4 ) ( 4 + j3 )

Fig. 2.28 Model with Divide Block

1.5 The Bus Creator and Bus Selector Blocks:

Bus Creator blocks to create signal buses and Bus Selector blocks to access the
components of a bus.

Example 2.5: Solve the following functions. The following model simulates the combined
functions of
a. sin2t
d
b. ( sin2t )
dt
c.  sin2t dt

Fig. 2.29 Model with Bus Creator

20
 Fig. 2.30 Output from the model with bus creator

The Bus Selector Model:

Double-click the bus selector block to select the output signal(s).

Fig. 2.31 Model with Bus Creator and Bus Selector Blocks

21
 Fig. 2.32 Model with Bus Selector Block

Fig. 2.33 Output from the model with bus selector block

The Mux and De-mux Blocks:

The Simulink Mux block combines its inputs into a single output. An input can be a scalar,
vector, or matrix signal. The Simulink Demux block extracts the components of an input signal
and outputs the components as separate signals. The block accepts either vector (1−D array)
signals or bus signals

Note: A sine wave can be converted into a Cosine wave by changing the phase angle to pi/2.

22
 Fig. 2.34 Model with Mux Block

Fig. 2.35 Output from the Mux Block

Double−click on the Constant Block and enter the given row vector [ 1 2 3 4 5 6 7 8 9]

Fig. 2.36 Constant block with (1x9) row vector

Double-click on the Demux block and for the Number of outputs enter [-1, 4, -1]

23
 Fig. 2.37 Model with Demux Block

Fig. 2.38 Output from Demux Block

The switch block:

The Switch block will output the first input or the third input depending on the value of the
second input. The first and third inputs are called data inputs. The second input is called the
control input.

Fig. 2.39 Switch Block set for u2>=0, Threshold = 0

24
 Fig. 2.40 Switch Block set for u2>=0, Threshold = 76

The Sum Block:

This block performs addition or subtraction on its inputs. This block can add or subtract scalar,
vector, or matrix inputs.

Example 2.6: Perform the operation A + B – C. Given,

A = [3+4j 1+0j 5-2j; 2-3j 4+j 7-4j; 1+6j 8-5j 4+7j];
B = [4+3j 0+2j -2+5j; -3+2j 6+7j -3-4j; 1+8j -5-3j 2-7j];
C = [-2+3j 7+2j -5-2j; 3-2j 4-7j -4+3j; -3+8j 7-4j -6+9j];

Fig. 2.41 Model for Example 2.6

25
The Gain Block:

The Gain block multiplies the input by a constant value (gain). The input and the gain can
each be a scalar, vector, or matrix. We specify the value of the gain in the Gain parameter.

Example 2.7:
Perform the matrix multiplication A  B given A = [1 -1 2] and B = [2 3 4]. Use an apostrophe
to convert A of 1x3 to 3x1. Double-click the gain block and enter the value of gain as k. In
MATLAB’s command window type k = [2 3 4].

Fig. 2.42 Model for Example 2.7

The Relational Operator Block:

The Relational Operator block performs the specified comparison of its two inputs. We select
the relational operator connecting the two inputs with the Relational Operator parameter.

Operation description:
== TRUE if the first input is equal to the second input
~= TRUE if the first input is not equal to the second input
< TRUE if the first input is less than the second input
<= TRUE if the first input is less than or equal to the second input
>= TRUE if the first input is greater than or equal to the second input
> TRUE if the first input is greater than the second input

Example 2.8: Determine the determinants of matrices A and B and check whether they are
equal or unequal.

26
  1 2 3 4 5 6
A =  4 5 6  B =  1 2 3 
7 8 9  7 8 9 

In MATLAB’s Command Window, enter A and B as follows:

>> A = [1 2 3;4 5 6;7 8 9];
>> B = [4 5 6;1 2 3;7 8 9];

Fig. 2.43 Model for Example 2.8

The Logical Operator Block:

The Logical Operator block performs the specified logical operation on its inputs. An input
value is TRUE (1) if it is nonzero and FALSE (0) if it is zero.

Table 1: Operation description

AND TRUE if all inputs are TRUE
OR TRUE if at least one input is TRUE
NAND TRUE if at least one input is FALSE
NOR TRUE when no inputs are TRUE
XOR TRUE if an odd number of inputs are TRUE
NOT TRUE if the input is FALSE and vice-versa


Example 2.9: Simulate the following Boolean expression D = A  B + C (where the dot

denotes the ANDing of the variables, A and B the plus (+) sign denotes the ORing of A  B
with C) using the Logical Operator Block. Let A = 1, B = 0, and C = 1.

27
 Fig. 2.44 Model for Example 2.9

The Saturation Block:

The Saturation block establishes upper and lower bounds for an input signal. When the input
signal is within the range specified by the Lower limit and Upper limit parameters, the input
signal passes through unchanged. When the input signal is outside these bounds, the signal
is clipped to the upper or lower bound.

Example 2.10: Perform the function y = 3x where x and y are specified in MATLAB’s
2

command window in the following format.

>> x = 0:1:6;
>> y = 3.*x.^2;

Double-click the saturation block to change the upper and lower limit value.

Fig. 2.45 Double-click on the Saturation Block

28
 Fig. 2.46 Model for Example 2.10 (Upper limit = 75 and Lower limit = 1)

The Integrator Block:

The Integrator block integrates its input and it is used with continuous−time signals.

➢ Double-click the integrator block to get the following window.

Fig. 2.47 Function Block Parameters for the Continuous Integrator Block

29
Table 2: Integrator Block Parameters

Integrator Block 1 Default (Initial condition source: internal)

Integrator Block 2 Initial condition source: external. All other parameters are in their
default states.
Integrator Block 3 Limit output: check mark. All other parameters are in their default
states.
Integrator Block 4 External reset: rising. All other parameters are in their default states.
Integrator Block 5 External reset: falling. All other parameters are in their default states.
Integrator Block 6 External reset: either. All other parameters are in their default states.
Integrator Block 7 External reset: level. All other parameters are in their default states.
Integrator Block 8 Show state port: check mark. All other parameters are in their default
states.
Integrator Block 9 External reset: rising.
Initial condition source: external.
Limit output: checkmark
Show saturation port: checkmark
Show state port: checkmark

Fig. 2.48 Different configurations for the Integrator Block

Example 2.11: Simulate the following differential equation subject to the initial conditions x(0)
= 0.5 and x’(0) = 0.

d2 x dx
2
+4 + 3x = 3  f ( t )
dt dt

The Constant 1 and Constant 2 blocks represent the initial conditions.

30
 Fig. 2.49 Model for Example 2.11

Algebraic Loop
An algebraic loop is formed when two or more blocks with direct feedthrough (the output of
the block at time t, is a function of the input at time t) form a feedback loop. The basic problem
with algebraic loops is that the output, y, at time, t, is a function of itself. An algebraic loop
generally occurs when an input port with direct feedthrough is driven by the output of the same
block, either directly, or by a feedback path through other blocks with direct feedthrough.

Example 12: Perform the Algebraic Loop y = 10 − y or y = 5

Fig. 2.50 Model for Example 2.12

The Integrator block's state port makes it possible to avoid creating algebraic loops when
creating an integrator that resets itself based on the value of its output.

The Unit Delay Block:

31
The Unit Delay block holds and delays its input by the sample period you specify. When placed
in an iterator subsystem, it holds and delays its input by one iteration. This block is equivalent
to the z-1 discrete-time operator. The block accepts one input and generates one output.

The Discrete–Time Integrator Block:

The Discrete−Time Integrator block performs discrete−time integration or accumulation of a

signal. We use this block in discrete−time systems instead of the Continuous Integrator block
in continuous−time systems. This block can integrate or accumulate using the Forward Euler,
Backward Euler, and Trapezoidal methods.

Example 2.13: Given the input: Pulse Type: Time-based; Time (t): Simulation time;
Amplitude: 0.25; Period (secs): 2; Pulse Width (% of Period): 50; Phase delay (secs): 0

Fig. 2.51 Model with Unit Delay and Discrete-Time Integrator Blocks

Output from Scope 1 Output from Scope 2

32
 Fig. 2.53 Output(s) from Example 2.13

Pulse Generator Block:

Fig. 2.54 Double-clicking the Pulse Generator Block

The period of the Pulse Generator Block depends on the switching frequency. If the switching
frequency is 50 Hz, then the period is (1/frequency) equal to 0.02. Pulse width (% of period)
represents the on and off time of the signal. In this case, the pulse width is taken as 50%.
Phase delay is entered as 0.01.

33
 Block Parameters: Amplitude = 1; Period = Block Parameters: Amplitude = 1; Period = 0.02
0.02 secs; Pulse width = 50 (% of Period); secs; Pulse width = 50 (% of Period); Phase
Phase delay = 0. delay = 0.01.

Fig. 2.55 Outputs from Pulse Generator Block with Different Block Parameters

The Vector Concatenate Block:

The Concatenate block concatenates the signals at its inputs to create an output signal whose
elements reside in contiguous locations in memory.

Example 2.14: Given A = [1 2 3 4] and B = [5 6 7 8]

Fig. 2.56 Model for Example 2.14

Simulink Library Browser: Simulink contains a large number of blocks from which models can be
built. These blocks are arranged in Block Libraries which are accessed in the Simulink library browser
window shown below. Each icon in the main Simulink window can be double-clicked to bring up the
corresponding block library. Blocks in each library can then be dragged into a model window to build a
model.

34
 Simulink

Continuous

Discontinuities

35
 Discrete

Logic & Bit

Operations

Math
Operations

36
 Sinks

Sources

Some Important Block(s):

From and Goto Blocks:

The From block accepts a signal from a corresponding Goto block

Fig. 2.57 Model Using From and Goto Blocks

37
From Workspace Block:

The From Workspace block reads data from the MATLAB workspace. The workspace data
are specified in the block's Data parameter via a MATLAB expression that evaluates to a 2−D
array.

Assignment 2
In this assignment, you will be required to complete the Simulink Onramp self-paced online
course offered by MathWorks. Simulink is a powerful simulation and modeling tool widely used
in engineering and scientific applications. This course will introduce you to the basics of
Simulink and guide you through various hands-on exercises to develop your skills. After
completing the course and obtaining the certificate, save it in a digital format (e.g., PDF or
image file). Name the file using the following format: [YourName_RegNo]. Submit your
certificate by uploading it to the assignment submission portal and a hard copy to the
instructor. The deadline for submitting the certificate is xx-August-2023.

38
 Tutorial 3: Solve the Systems of Linear Equations

Using Algebraic Constraint blocks found in the Math Operations library, Display blocks
found in the Sinks library, and Gain blocks found in the Commonly Used Blocks library,
create a model that will produce the simultaneous solution of three equations with three
unknowns.
The model will display the values for the unknowns z1 , z 2 and z 3 in the system of the

equations
a1z1 + a2 z 2 + a3 z3 + k1 = 0
a 4 z1 + a5 z 2 + a6 z3 + k 2 = 0
a7 z1 + a8 z 2 + a9 z3 + k 3 = 0
Simulink Model:

Fig. 3.1 Model for Tutorial 3

Go to MATLAB’s Command Window and enter the following values:

a1 = 2; a2 = −3; a3 = −1; a4 = 1; a5 = 5; a6 = 4; a7 = −6; a8 = 1; a9 = 2

k1 = −8; k 2 = −7; k 3 = 5

39
 Tutorial 4: Modeling and Simulation of Mass-Spring-Damper System
Using Simulink

Example 4.1: Simulate the Mass-Spring-Damper System using Gain Blocks and Integrators.

Fig. 4.1 Free-Body Diagram of Mass-Spring-Damper System

Equation of Motion:

d2 x dx
m +c + kx = F ( t )  m x + c x + kx = F ( t )
dt dt

Inputs: Mass (m) = 1 kg; Damping constant (c) = 0.2 Ns/m; Spring constant (k) = 1 N/m
Force Input (F) = 1 N (It is a step function with magnitude 1 at t = 0)

• Set the block parameters of the step input

➔ Step time = 0
➔ Initial value = 0
➔ Final value = 1
• Initial conditions: x(0) = 0 and x’(0) = 0.

Start by solving the model equation for the highest-order derivative term.

d2 x dx
m 2 = F(t) − c − kx
dt dt

The first block you create should be a Sum block, where the output of the Sum block is the
left-hand term of the equation above. The Sum block is to have 3 inputs. Double-click on the
Block and set the parameters to rectangular and + - -.

Fig. 4.2 Sum Block with Output Line

40
Add a gain (multiplier) block to normalize the coefficient, m, to modify the signal so it is equal
to the highest-order derivative term alone.

Fig. 4.3 Sum Block with Gain Block

Add two integrators to the model to obtain the velocity (dx/dt) and displacement (x).

Fig. 4.4 Sum Block with Gain and Integrator Blocks

In the model, connect the integrated signals with gain blocks to create the terms on the right-
hand side of the model equation.

Fig. 4.5 Sum Block with Gain and Integrator Blocks

Connect all the input signals to the appropriate inputs of the Sum Block.

41
 Fig. 4.6 Model Representing Mass-Spring-Damper System (Closed Loop)

Add a Step Block for force input from the Source library and a Scope Block from the Sink
library.

Fig. 4.7 Model Representing Mass-Spring-Damper System without Block Parameters

Final Simulink model for a mass-spring-damper system with input values.

Fig. 4.8 Model Representing Mass-Spring-Damper System with Block Parameters

42
To pass the output to MATLAB’s workspace, add them To Workspace block to the line as
shown in Figure.

Fig. 4.9 Model for Mass-Spring-Damper System with Output Line to Workspace

Fig. 4.10 Output from the MSD Model

Fig. 4.11 MSD Model without Block Parameters

43
MATLAB Code to Run the Model Shown in Figure 4.11:

➢ Create a new script file in MATLAB to write this code.

clc;
clear;
close all;

%% Input parameters

m = 1; % mass (kg)
c = 0.2; % damping coefficient (Ns/m)
k = 1; % spring constant (N/m)
TF = 50; % Stop time (sec)

%% To run the Simulink model

Result = sim("Mass_Spring_Damper_Model.slx");

%% To plot the output

plot(Result.simout,'-r','linewidth',2)
title('Mass-Spring-Damper System')
ylim([0 2.0])
ylabel('Displacement (m)')

grid on
ax = gca;
ax.GridColor = [0 0 0];
ax.GridLineStyle = ':';
ax.GridAlpha = 0.3;
ax.Layer = 'bottom';

%% To save the figure in png format

print('MSD','-dpng','-r600');

44
Example 4.2: Mass-Spring-Damper (MSD) System Modeling Using State-space Model and
Transfer Function Model in Simulink

State-space Model:
To determine the state-space representation of the mass-spring-damper system, we must
reduce the second-order governing equation to a set of two first-order differential equations.
To this end, we choose the position and velocity as our state variables.

 X   x  Displacement 
X =  1 =   =  
 X2   x   Velocity 

From the Equation of Motion for Mass-Spring-Damper System:

d2 x 1  1
− kx  x =  F ( t ) − c x − kx  = (F ( t ) − cX2 − kX1 )
dx
m = F(t) − c
dt 2
dt m  m

The state equation in this case is:

     X2 
 X1   x   Velocity   
X= = = =
     Acceleration  1 (F ( t ) − cX2 − kX1 ) 
 X2   x  m 

 X2 
= 0 1   X1   0 
X=1 + F(t)
 (F ( t ) − cX2 − kX1 )   −k m −c m  X2  1 m

m 

If, for instance, we are interested in controlling the position of the mass, then the output
equation is:

X 
y = 1 0  1 
 X2 

MATLAB Code:

m = 1;
k = 1;
c = 0.2;
F = 1;
A = [0 1; -k/m -b/m];
B = [0 1/m]';
C = [1 0];
D = [0];
sys = ss(A,B,C,D)

45
The State-Space Block:

The State-Space block implements a system defined by the state−space equations

x = Ax + Bu

y = Cx + Du

where x and u are column vectors, matrix A must be an n x n square matrix where n represents
the number of the states, matrix B must have dimension n x m where m represents the number
of inputs, matrix C must have dimension r x n where r represents the number of outputs, and
matrix D must have dimension r x m.

A = [0 1; -k/m -c/m];
B = [0 1/m]';
C = [1 0];
D = [0];

Double-click the state-space model and set the block parameters as shown in Figure 4.12.

Fig. 4.12 State-Space Block Parameters

46
 Fig. 4.13 Model with State-Space Model

• Set the block parameters of the step input

➔ Step time = 0
➔ Initial value = 0
➔ Final value = 1

Fig. 4.14 Output from State-Space Model

The Transfer Function Block:

The Transfer Fcn block implements a transfer function where the input F(s) and output H(s)
can be expressed in transfer function form as the following expression.

Output H ( s )
G (s) = =
Input F (s)

Transfer Function Representation:

we can determine the transfer function directly from the state-space representation as follows:

47
 H(s)
G ( s) = = C ( sI − A ) B + D
−1

F ( s)

The Laplace transform for this system assuming zero initial conditions is

ms2 X ( s ) + csX ( s ) + kX ( s ) = F ( s )

and, therefore, the transfer function from force input to displacement output is

X ( s) 1
=
F (s) ms + cs + k
2

MATLAB Code:

s = tf('s');
sys = 1/(m*s^2+b*s+k);

Alternate MATLAB Code:

num = [1];
den = [m b k];
sys = tf(num,den);

Fig. 4.15 Simulink Model for MSD System with Subsystem Block, State-Space Block,
and Transfer Function Block

48
 Fig. 4.16 Model Inside Mass_Spring_Damper Subsystem in the above model

Fig. 4.17 Block Parameters for State-Space Block (left) and Transfer FCN Block

Output from Scope 1, 2, and 4 Output from Scope 3

49
Example 4.3: Modeling of MSD System Using MATLAB Function Block

The MATLAB FCN Block:

The MATLAB Fcn block applies the specified MATLAB function or expression to the input.

Fig. 4.19 MSD System Using MATLAB FCN Block (Inside MSD System)

Fig. 4.20 Model Inside MSD System in Fig. 4.19

MATLAB Code (Inside MATLAB FCN Block)

function y = fcn(F,dx_dt,x)
% m = 100; b = 300; k = 6000
y = (F/100)-((300*dx_dt)/100)-((6000*x)/100);

50
Example 4.4: Mass-Spring-Damper Model Using Integrator, Second-Order Block

The Integrator, Second-Order Block:

The Second-Order Integrator block and the Second-Order Integrator Limited block solve the
second-order initial value problem:

d2 x
= u; ( dx dt )t =0 = dx 0 ; ( x )t =0 = x 0
dt 2

SIMULATION > Library Browser > Simulink > Continuous > Integrator, Second-Order
Block

Fig. 4.21 Model for Example 4.4

Fig. 4.22 Output from the above Model

51
Example 4.5: Modeling of MSD System Using Fcn Block

Fcn Block:

The Fcn block applies a specified expression to its input denoted as u. If u is a vector, u(i)
represents the ith element of the vector; u(1) or u alone represents the first element. The
specified expression can consist of numeric constants, arithmetic operators, relational
operators, logical operators, and math functions.

Fig. 4.23 MSD Model with Fcn Block

Input Values for MSD System: m = 1 kg; c = 0.2 Ns/m; k = 1 N/m; F = 1 N

 
Function: x = (1 m )  F − c x − kx 
 

Fig. 4.24 Fcn Block Parameters

52
 Assignment 3

1. Simulate the following equations

M1a1 = F − k ( x1 − x 2 ) − μM1g x1

M2a2 = k ( x1 − x 2 ) − μM2 g x 2

Given: Mass_1 = 1 kg; Mass_2 = 0.5 kg; k = 1 N/sec; F = 1 N; mu = 0.02 sec/m; g = 9.8 m/s2.
Use Signal Generator Block to input the force (F). set Amplitude = -1 and Frequency = 2 Hz.
Set the stop time of the simulation to 1000 sec.

The Signal Generator Block:

The Signal Generator block can produce one of four different waveforms: sine wave, square
wave, sawtooth wave, and random wave. The signal parameters can be expressed in Hertz
(the default) or radians per second. We can invert the waveform by specifying a negative
amplitude in the block’s parameters window.

Set the frequency to 0.2 Hz for all the Signal Generator Blocks.

53
 2. Simulate the following model using Simulink:

Equations of Motion:

m1a1 = −b1(x1 − x 2 ) − k1 ( x1 − x 2 ) + u
 
m2a2 = b1(x1 − x 2 ) + k1  w − x 2  + k 2 ( w − x 2 ) − u
 
u and w are control inputs. Take step input for u (Step Time to "0" and leave its Final Value
"1") and w (Step Time to "0" and Final Value to "0").

The system parameters are as follows:

(m1) Body mass 2500 kg
(m2) Suspension mass 320 kg
(k1) Spring constant of suspension system 800000 N/m
(k2) Spring constant of wheel and tire 500000 N/m
(b1) Damping constant of suspension system 350 N.s/m
(b2) Damping constant of wheel and tire 15020 N.s/m

Fig. 4.28 Output of the Suspension Model

54
 Tutorial 5: Interaction of a Simulink Model from a MATLAB Script

Example 5.1: This example illustrates how to control and interact with a Simulink model from
a MATLAB script. This is useful if you would like to analyze data generated from a Simulink
model in the MATLAB environment.

Simple Pendulum Model:

Equation of Motion:
 1  g  c 
θ(t) =  2 
T(t) −   θ ( t ) −  2  ( )
θ t
 mL  L  mL 

where T = Torque; c = Damping constant

Model Input Parameters:

Mass = m = 1 kg
Acceleration due to gravity = g = 9.81 m/s2
Damping constant = c = 0.2 Ns/m
Length = L = 1m

Fig. 5.1 Closed Loop Model of Simple Pendulum

55
Defining Simulink model parameters in the MATLAB workspace:

clear
clc
close all

%% Define Model Constants

g = 9.81; % Acceleration due to gravity

m = 1; % Mass of the bob
c = 0.2; % Damping constant
L = 1; % Length of the rod

Fig. 5.2 Model with Step Input and Scope Blocks

%% Step Input Time and Stop Time Setting

Step_Time = 10; % step Time in seconds

tFinal = 50; % Stop Time in seconds

Fig. 5.3 Step (Torque) Input Block Parameters (Double-click the Step Block)

56
SIMULATION > Stop Time (enter → tFinal); Default value = 10 seconds

Fig. 5.4 Stop Time Setting

%% Initial Condition for the Integrators

thetaDot0 = 0;
theta0 = 0;

Fig. 5.5 Block Parameters: Integrator Block

Running a Simulink model from a MATLAB script:

➢ Save the Simulink file as Simple_Pendulum (file will be saved as

Simple_Pendulum.slx)
➢ Save the MATLAB script file as Run_Simple_Pendulum

Note: Both files should be in the same folder.

%% To Run the Simulink File

Result = sim("Simple_Pendulum.slx");

57
Sending simulation data back to the MATLAB workspace:
Using Out1 Block:

To activate the signal dimension:

Display the dimensions directly on the block diagram. Use this technique to trace signal
dimensions along a path of blocks. In the model, on the Debug tab, select Information
Overlays > Signal Dimensions).

%% Extract the Data Generated Using Out1 Block

thetaSignal = Result.yout.getElement('Out');
time = Result.tout;
thetaDot = thetaSignal.Values.Data(:,1);
theta = thetaSignal.Values.Data(:,2);

To Find the Peak Value of Theta and the Corresponding Time:

[thetaMax,indexMax] = max(theta);
tMax = time(indexMax);

Plot a Graph Using Generated Data:

figure;

subplot(2,1,1)
hold on
plot(time,theta,'r-','linewidth',1.5);
plot(tMax,thetaMax,'bx','linewidth',2,...
'markersize',14)
grid on;
xlabel('Time (seconds)');
ylabel('Theta (rad)');
legend('Theta','MaxTheta')
title('Data Generated by the Simulink Model')
subplot(2,1,2)
plot(time,thetaDot,'m:','linewidth',1.5)

58
grid on;
xlabel('Time (seconds)');
ylabel('ThetaDot (rad/s)');

%% To save the figure in png format

print('Simple_Pendulum','-dpng','-r600');

To make the curve smooth in the above plot: SIMULATION > PREPARE > Model Settings
(drop down) > Solver > Max step size (enter Max_Step) > Ok

59
Add a line in MATLAB Script (before Result = sim("Simple_Pendulum.slx"))

%% Maximum Step Size

Max_Step = 0.01;

Using MATLAB data as input to a Simulink model:

Add a line in MATLAB Script (before Result = sim("Simple_Pendulum.slx"))

%% From Workspace

xT = [t u];

Note: Create the column vectors for t and u based on the requirement.

60
 Tutorial 6: Modeling and Simulation of a DC Motor Using Simulink

This tutorial illustrates the modeling of a DC motor by summing the torques acting on the rotor
inertia and integrating the acceleration to give velocity.

The physical parameters for this example are:

(J) Moment of inertia of the rotor 0.01 kgm2

(b) Motor viscous friction constant 0.1 N.m.s
(Ke) Electromotive force constant 0.01 V/rad/sec

(Kt) Motor torque constant 0.01 N.m/Amp

(R) Electric resistance 1 Ohm

(L) Electric inductance 0.5 H
Applying Newton's law and Kirchoff's law to the motor system to generate the following
equations:

d2θ dθ d2θ 1  dθ 
J 2 = T −b  2 =  K ti − b 
dt dt dt J dt 

di di 1  dθ 
L = −Ri + V − e  =  V − Ke − Ri 
dt dt L  dt 

Simulink Model:

Fig. 6.1 Model for Tutorial 6

61
 Tutorial 7: Design and Implementation of PID Controller in Simulink

The PID Controller block implements a PID controller (PID, PI, PD, P only, or I only). The block
is identical to the Discrete PID Controller block with the Time domain parameter set
to Continuous-time.

We can also model the PID controller using the blocks shown in the figure below.

Fig. 7.1 Simulink Model for PID Controller

Example 7.1: Simulate a second-order differential equation from Mass-Spring-Damper

System with a PID controller. The equation model is given directly as a transfer function.
Develop a model using the PID controller block. Use the Signal Generator Block as a forcing
input with Waveform of a square signal, Time(t) of Use simulation time, Amplitude of 1, and
Frequency of 1 rad/sec. Also, P – I – D tuning parameters are also shared to avoid spending
time on the tuning of controllers (Kp = 300; Ki = 350; Kd = 10).

H( s) 20
G ( s) = =
F (s) s + 10s + 20
2

62
 Fig. 7.2 Model using the PID Controller Block

Example 7.2: Develop the same model given in Example 7.1 using individual gain blocks for
P – I – D.

Fig. 7.3 Implementing the Control Feedback System of a PID Controller Using Gain
Blocks

Fig. 7.4 Output from the Model Shown in Fig. 7.3

63
Example 7.3: Simulate a second-order differential equation from Mass-Spring-Damper
System with a PID controller. Use the Scope Block and dashboard scope block to visualize
the result.
H( s) 1
G ( s) = =
F (s) s + 5s + 6
2

Dashboard Scope Block:

Blocks that can control parameter values and display signal values during simulation.

Double-clicking the block, a window pop-up. Select the actual output line and checkbox to set
the output signal and then select the controlled output line and checkbox as shown in Fig. 00.

64
 Fig. 7.7 Simulation Settings

Fig. 7.8 Simulation Settings

Double-click the Scope Block to open the following window. You can run the simulation
directly from the Scope Block.

Fig. 7.9 Output of the Model for Example 7.3

To analyze the signal, you can use the cursor measurement option in the Scope.

65
 Fig. 7.10 Scope Block (On Double-clicking)

Double-click the PID Controller Block and Click the Tuner to tune the P, I, and D value.

Fig. 7.11 PID Controller Block Parameter

66
Adjust the Response Time slider and Transient Behavior slider to tune the PID value.

Fig. 7.12 PID Tuner Window

After achieving the desired response, click the Update Block to update the P, I, and D
value of the controller.

Fig. 7.13 PID Tuner

Fig. 7.14 PID Tuner

67
The data can be inspected using Data Inspector when using Dashboard Scope for
visualizing the data as shown in Fig. 7.15.

Fig. 7.15 Data Inspector

Example 7.4: Aircraft Pitch System Modeling

Simulate a second-order differential equation from aircraft pitching motion with a PID
controller. An equation model is given along with a state-space model and transfer function.
Develop a model using the PID controller block. Tune the PID controller to meet the following
design requirements:

• Overshoot less than 10%

• Rise time less than 2 seconds
• Settling time less than 10 seconds
• Steady-state error of less than 2%

Equations of Motion:

α = −0.313α + 56.7q + 0.232δ

q = −0.0139α − 0.426q + 0.0203δ

θ = 56.7q

These values are taken from the data from one of Boeing's commercial aircraft.

Transfer Function:
θ ( s) 1.151s + 0.1774
P (s) = =
Δ (s) s + 0.739s2 + 0.921s
3

State-Space Representation:

 
α   −0.313 56.7 0   α   0.232 
      
 q  =  −0.0139 −0.426 0   q  + 0.0203  δ
   56.7 0   θ   0 
θ  
0
 

68
 α 
y = 0 0 1  q 
θ 
MATLAB Code:

s = tf('s');
P_pitch = (1.151*s+0.1774)/(s^3+0.739*s^2+0.921*s);
t = [0:0.01:10];
step(0.2*P_pitch,t);
axis([0 10 0 0.8]);
ylabel('pitch angle (rad)');
title('Open-loop Step Response');

Open Loop System:

Fig. 7.16 Open Loop Model

Closed Loop System:

Fig. 7.17 Closed Loop Model

Closed Loop System with Gain constant:

Fig. 7.18 Closed Loop Model with Gain Constant

69
 Fig. 7.20 Plant Model with Gain as Controller

The gain was designed using the Linear Quadratic Regulator method and resulted in a

calculation of K = [-0.6435 169.6950 7.0711].

Output from Open Loop Model Output from Controlled Model

Fig. 7.22 Model with Pre-compensator Constant (Nbar = 7.0711)

The steady-state error requirement is not met since the response does not settle to within 2%
of the commanded reference of 0.2 radians. This deficiency was addressed by adding a
constant pre-compensator = 7.0711 to scale the output to the desired level.

70
 Fig. 7.23 Output from the Model Shown in Fig. 7.22

In order to investigate the efficiency of the pre-compensator, let's add a disturbance to our
model as shown in the figure below. The disturbance is generated by a Step block with
the Final value set to "0.2" and the Step time set to "3".

Fig. 7.24 Model with Added Disturbance

Fig. 7.25 Output from the Model Shown in Fig. 7.24

71
 Fig. 7.26 Model with PID Controller

➢ Double-click the PID Controller Block

Fig. 7.27 PID Controller Block Parameter

Adjust the Response Time (seconds) and Transient Behavior to improve the performance of
the system.

Fig. 7.28 PID Tuner

Click Show Parameters shown in Fig. 7.28 to see the final parameters of the system
shown in Fig. 7.30.

72
Response of the System after tuning PID values.

Fig. 7.29 Output After PID Tuning

Fig. 7.30 Performance Parameters

Example 7.5: Design PID Controller for Aircraft Pitch System Modeling

Kp = 10.5244070119708; Ki = 4.4379063209003; Kd = 5.55492229740641;

Filter coefficient = 760.080531029323

73
 Fig. 7.31 Model with Designed PID Controller

Fig. 7.32 Model with In-built PID Block

Fig. 7.33 Output from Model with PID Controller

Example 7.6: Extract the given model into MATLAB (linearization)

Fig. 7.34 Closed Loop Model

The Simulink Control Design toolbox offers the functionality to extract a model from Simulink
into the MATLAB workspace.

74
Right-click on the step input signal (output of the Step block) and choose Linear Analysis
Points > Input Perturbation from the resulting menu to identify the input of our closed-loop
system.

Next, right-click on the output line of the transfer function and select Linear Analysis Points
> Output Measurement from the menu to choose the output of our system.

Your model should now appear as follows where the small arrow symbols identify the input
and output of the model.

We can now extract the model by opening the Linear Analysis Tool. This is accomplished
by selecting Model Linearizer from under the Apps menu at the top of the model window.
Following these steps will open the window shown below.

This tool generates an LTI object from a (possibly nonlinear) Simulink model and allows you
to specify the point about which the linearization is performed. Since our Simulink model is
already linear, our choice of operating point will have no effect and we can leave it as the
default Model Initial Condition. In order to generate the linearized model, select
the Step button in the above figure, which is indicated by a small green triangle. The Linear
Analysis Tool window should now appear as shown below.

75
 Right-click on the linsys1 under Linear Analysis Workshop to export the
linear model to MATLAB Workspace

On right-clicking

Example 7.7: Simulate and test a controller for a solar panel as it tracks the movement of the
sun throughout the day. Take the stop time for this problem as 15 hours.

Equation of Motion for the Solar Panel:

d2θ 1  dθ 
=  T − Kd
dt 2
J dt 
where T = Torque; dθ/dt = Angular velocity; J = Rotational Inertia = 0.01 kg/m2; Kd = Constant
= 0.1.

76
Equation of Motion for the DC Motor:

di 1  dθ 
=  V − K gK f − Ri 
dt L  dt 

T = K gK ti

where L = Electric Inductance = 0.5 H; V = Voltage; Kg, Kf and Kt = Constant = 0.01; i =

Armature current; Electric resistance = 1 Ohm

To generate Sun position data in the MATLAB workspace:

MATLAB Code:

a1 = 422.7;
b1 = 4.954e-05;
c1 = 0.6253;
a2 = 420.5;
b2 = 4.976e-05;
c2 = 3.768;
a3 = 0.1964;
b3 = 0.0002443;
c3 = -0.08255;
x = (0:360:54000);
f = (a1*sin(b1*x+c1) + a2*sin(b2*x+c2) + a3*sin(b3*x+c3));
sunTime = x';
sunPosition = f';

Fig. 7.40 Model for Solar Panel

Fig. 7.41 Model for DC Motor

77
 Fig. 7.42 Model without Controller

Fig. 7.43 Output from Model without Controller

Fig. 7.44 Model with PID Controller

78
Subsystem:

Fig. 7.45 Model Inside Subsystem

To import the Sun Time and Sun Position data:

SIMULATION > PREPARE > Click the drop-down > CONFIGURATION & SIMULATION >
Model Settings

By Clicking the Model Settings, a window pop-up as shown in the figure below:
➢ Select the Input and enter [sunTime sunPosition]
Note: The sunTime and sunPosition data must be in the workspace of MATLAB while
running this model.

79
Fig. 7.47 Output from the Model Shown in Fig. 7.44

80
 Tutorial 8: Extraction of a Linear Model into MATLAB

A linear model of the system can be extracted from the Simulink model into the MATLAB
workspace. This can be accomplished by employing the MATLAB command linmod or from
directly within Simulink.

Example 8.1: Develop a non-linear 6 DOF aircraft model using RCAM Model.

Description of Non-linear 6DOF Aircraft Model:

The aircraft is a twin-engine, commercial transport aircraft (Similar to Boeing 757-200).

Aircraft Specifications:
Aircraft total mass = 120000 kg; Mean aerodynamic chord (c_bar) = 6.6 m; Tail moment arm
(l_t) = 24.8 m; Wing planform area (S) = 260 m2; Tail planform area (St) = 64 m2;

Control Vectors:

 u1   δa   Aileron 
     
 u2   δe   Elevator 
u =  u3  =  δr  =  Rudder 
     
 u4   δth1   Throttle1 
 u   δ   Throttle2 
 5   th2   

State Vectors:

 x1   u   Velocity in body x − axis 

     
 x2   v   Velocity in body y − axis 
 x3   w   Velocity in body z − axis 
     
 x 4   p   Angular rate about body x − axis 
x =  x 5  =  q  =  Angular rate about body y − axis 
     
 x 6   r   Angular rate about body z − axis 
     
 x7   Φ   Bank Euler Angle

x   θ   Pitch Euler Angle 
 8    
 x  Ψ  Azimuth Euler Angle 
 9

81
Required Models:
1. Control limit/saturation
2. Intermediate variables
3. Non-dimensional aerodynamic force coefficients in Fs
4. Aerodynamic force in Fb
5. Non-dimensional aerodynamic moment coefficient about aerodynamic center in Fb
6. Aerodynamic moment about the aerodynamic center in Fb
7. Aerodynamic moment about the center of gravity in Fb
8. Propulsion
9. Gravity effects
10. Explicit First Order Form
Control limit/saturation:

Rate limit for throttle movement: -1.6° to 1.6°

Saturation limit for throttle movement: 0.5° to 10°

Rate limit for aileron deflection: -25° to 25°

Saturation limit for aileron deflection: -25° to 25°

Rate limit for tail-plane deflection: -15° to 15°

Saturation limit for tail-plane deflection: -25° to 10°

Rate limit for rudder deflection: -25° to 25°

Saturation limit for rudder deflection: -30° to 30°

Intermediate Variables:

➢ Flight speed = Vair = x12 + x 22 + x 32

82
 ➢ Angle of Attack = α = a tan2 ( x3 ,x1 )
➢ β = sin−1 ( x 2 Vair ) = Sideslip angle
➢ Q = 0.5  ρ  Vair
2
= Dynamic Pr essure
➢ ω = p q r  = Aircraft Angular Rates
T

V body = u v w  =  x1 x 2 x 3  = Total Velocity of Aircraft

T T
➢

Non-dimensional aerodynamic force coefficients in Fs:

 n ( α − αL =0 ) if α  14.5  ( π 180 )
CL,wb = 
a3 α + a2α + a1α + a0 Otherwise
3 2

n = 5.5; αL=0 = −11.5  ( π 180 ) ; a0 = 15.212; a1 = −155.2; a2 = 609.2; a3 = −768.5

Non-dimensional aerodynamic force coefficients in Fs:

The angle of attack of the tail is expressed as,

αt = α − ε + u2 + 1.3  x 5  (lt Vair )

ε = ( ε α )  α − αL =0  = 0.25 α + 11.5  ( π 180 ) 

u2 = δe = Elevator deflection
1.3  x 5  ( lt Vair ) = 1.3  q  (lt Vair ) = Dynamic Pitch Re sponse
CL,tail = 3.1 ( St S )  αt

Non-dimensional aerodynamic force coefficients in Fs:

Total drag coefficient

CD = 0.13 + 0.07 ( 5.5α + 0.654 )

2

Total side-force coefficient

CY = −1.6β + 0.24u3
 −CD  0.13 + 0.07 ( 5.5α + 0.654 ) 
2

 
CF =  CY  = 
S
−1.6β + 0.24u3 
 −CL   CL,wb + CL,tail 
 

 cosβ sinβ 0 
C = T ( β ) C =  − sinβ cosβ 0  CF
W S S
F F

 0 0 1

Aerodynamic force in Fb:

FA =  −D; Y; − L =  −CDQS; CYQS; − CLQS

S

83
 cos α 0 − sinα  −Dcos α + L sinα 
b 
F = T  α F =  0
S
1 
0  FA = 
S
FY 
A A

 sinα 0 cos α   −Dsinα − L cos α 

Non-dimensional Aerodynamic moment coefficient in Fb:

 −1.4β   −11 0 5 
 Cl,ac     
  Sl   S l2 
=  −0.59 − 3.1 t t  ( α − ε )  +  0 −4.03  t t  0 
b C
Cm,ac = Cm,ac  +
  SC     SC   V
 Cn,ac      air
body
 1 − α (180 15π ) β   1.7 0 −11.5 
 −0.6 0 0.22   u 
   1
 0 −3.1 S tlt SC ( )
0  u2 
 
 0 0 −0.63  u3 

Aerodynamic moment about C.G. in Fb:

b b
Mac = Cm,ac QSC

( )
b b
Mcg = Mac + F Aero  r cg − r ac
T
r cg =  Xcg Ycg Zcg  = 0.23C 0 0.1C 

r ac =  Xac Zac  = 0.12C 0 0 

T
Yac

Propulsion Effects:

Fi = δth,i  m  g
F1 = u4  m  g
F2 = u5  m  g
At max imum thrust (u4 = u5 = 10 ( π 180 ) )
F1,max + F2,max T
= = 0.35
mg W

FE,i = Fi 0 0 body

b T

b b b
FE = FE,1 + FE,2
b b b
ME,cg = ME,i  FE,i
b T
ME,i =  Xcg − X APT,i YAPT,i − Ycg Zcg − Z APT,i 
b b b
ME,cg = ME,1 + ME,2

XAPT,1 YAPT,1 Z APT,1 XAPT,2 YAPT,2 ZAPT,2

Engine
Parameter
0.0 m -7.94 m -1.9 m 0.0 m 7.94 m -1.9 m

84
Gravity Effects:

0  0 
FG,E 
=0  =  0 
 W Earth mgEarth
1 0 0  cosΘ 0 − sinΘ   cosΨ sinΨ 0   0 
FG,B 
= 0 cosΦ sinΦ   0 1 0   − sinψ cosΨ 0   0 
0 − sinΦ cosΦ   sinΘ 0 cosΘ   0 0 1 mgE
 −mgsin ( x 8 ) 
 
FG,B =  mgsin ( x 7 ) cos ( x 8 ) 
mgcos ( x 7 ) cos ( x 8 ) 
 B
Explicit First Order Form:

   
x  u 
    1 body body body
y =  v  = F −ω V
    m
z w 
   
 40.07 0 2.098 
Ib =  0 64 0 
 2.098 0 99.92 
body body body body
F = Fgravity + Fengine + Faero
   
 x 4  p 
    −1 b
(b
 x 5  = q = Ib Mcg − ω  Ib ω
   
b
)
x6   r 
   
b body body
Mcg = Mengine + Maero
Explicit First Order Form:

   
 x 7  Φ  1 tanΘ sinΦ tanΘcosΦ  p 
    
 x 8  =  Θ  = 0 cosΦ − sinΦ  q
    0 sec Θ sinΦ sec ΘcosΦ   r 
 x 9  Ψ   
   
Blocks Required:

• Interpreted MATLAB Function

• Integrator
• Constant

85
 Fig. 8.1 Non-Linear 6 DOF Aircraft Model (RCAM Model)

MATLAB Code:

function[XDOT]=Nonlinear_6DOF_Model(X,U)

% State vector

x1 = X(1);
x2 = X(2);
x3 = X(3);
x4 = X(4);
x5 = X(5);
x6 = X(6);
x7 = X(7);
x8 = X(8);
x9 = X(9);

u1 = U(1);
u2 = U(2);
u3 = U(3);
u4 = U(4);
u5 = U(5);

%% Constants

m = 120000;
cbar = 6.6;
lt = 24.8;
S = 260;
St = 64;
Xcg = 0.23*cbar;
Ycg = 0;
Zcg = 0.10*cbar;

Xac = 0.12*cbar;
Yac = 0;
Zac = 0;

86
Xapt1 = 0;
Yapt1 = -7.94;
Zapt1 = -1.9;

Xapt2 = 0;
Yapt2 = 7.94;
Zapt2 = -1.9;

%% other constants

rho = 1.225;
g = 9.81;
depsda = 0.25;
alpha_L0 = -11.5*(pi/180);
n = 5.5;
a3 = -768.5;
a2 = 609.2;
a1 = -155.2;
a0 = 15.212;
alpha_switch = 14.5*(pi/180);

%---------------------Intermediate Variables-----------------%

% Calculated airspeed

Va = sqrt(x1^2+x2^2+x3^2);

% Calculate alpha and beta

alpha = atan2(x3,x1);
beta = asin(x2/Va);

% Calculate dynamic pressure

Q = 0.5*rho*Va^2;

% Also define the vector omega and v_b

omega_b = [x4;x5;x6];
V_b = [x1;x2;x3];

%------------------Aerodynamic Force Coefficient-------------%

if alpha <= alpha_switch

CL_wb = n*(alpha-alpha_L0);
else
CL_wb = a3*alpha^3+a2*alpha^2+a1*alpha+a0;
end

% Calculate lift coefficient for tail

87
epsilon = depsda*(alpha-alpha_L0);
alpha_t = alpha-epsilon+u2+1.3*x5*(lt/Va);
CL_t = 3.1*(St/S)*alpha_t;

% Total lift force

CL = CL_wb + CL_t;

% Total Drag Force

CD = 0.13+0.07*(5.5*alpha+0.654)^2;

% Calculate side force

CY = -1.6*beta + 0.24*u3;

%-----------------Dimensional Aerodynamic Forces-------------%

FA_s = [-CD*Q*S;
CY*Q*S;
-CL*Q*S];

% Rotate the force to body axis

C_bs = [cos(alpha) 0 -sin(alpha);

0 1 0;
sin(alpha) 0 cos(alpha)];

FA_b = C_bs*FA_s;

%--------Aerodynamic moment about Aerodynamic center---------%

eta11 = -1.4*beta;
eta21 = -0.59-(3.1*(St*lt)/(S*cbar))*(alpha-epsilon);
eta31 = (1-alpha*(180/(15*pi)))*beta;

eta = [eta11;
eta21;
eta31];

dCMdx = (cbar/Va)*[-11 0 5;
0 (-4.03*(St*lt^2)/(S*cbar^2)) 0;
1.7 0 -11.5];

dCMdu = [-0.6 0 0.22;

0 (-3.1*(St*lt)/(S*cbar)) 0;
0 0 -0.63];

% Calculate CM about aerodynamic center in Fb

88
CMac_b = eta + dCMdx*omega_b + dCMdu*[u1;u2;u3];

%----------------Aerodynamic Moment about AC-----------------%

% Normalized

MAac_b = CMac_b*Q*S*cbar;

%----------------Aerodynamic Moment about CG-----------------%

rcg_b = [Xcg;Ycg;Zcg];
rac_b = [Xac;Yac;Zac];
MAcg_b = MAac_b + cross(FA_b,rcg_b - rac_b);

%-----------------Engine Force and Moment--------------------%

F1 = u4*m*g;
F2 = u5*m*g;

% Assume that the engine thrust is aligned with Fb

FE1_b = [F1;0;0];
FE2_b = [F2;0;0];

FE_b = FE1_b + FE2_b;

% Engine moment due to offset of engine thrust from CoG

mew1 = [Xcg - Xapt1;

Yapt1 - Ycg;
Zcg - Zapt1];

mew2 = [Xcg - Xapt2;

Yapt2 - Ycg;
Zcg - Zapt2];

MEcg1_b = cross(mew1,FE1_b);
MEcg2_b = cross(mew2,FE2_b);

MEcg_b = MEcg1_b + MEcg2_b;

%--------------------Gravity Effects-------------------------%

g_b = [-g*sin(x8);
g*cos(x8)*sin(x7);
g*cos(x8)*cos(x7)];

Fg_b = m*g_b;

%------------------------State Derivative--------------------%

89
% Inertia Matrix

Ib = m*[40.07 0 -2.0923;
0 64 0;
-2.0923 0 99.92];

% Inverse of Inertia Matrix

invIb = (1/m)*[0.0249836 0 0.000523151;

0 0.015625 0;
0.000523151 0 0.010019];

% All the forces in Fb and calculate udot, vdot, wdot

F_b = Fg_b + FE_b + FA_b;

x1to3dot = (1/m)*F_b - cross(omega_b,V_b);

% All the moments in Fb and calculate pdot, qdot, rdot

Mcg_b = MAcg_b + MEcg_b;

x4to6dot = invIb*(Mcg_b - cross(omega_b,Ib*omega_b));

% Calculate phidot, thetadot, psidot

H_phi = [1 sin(x7)*tan(x8) cos(x7)*tan(x8);

0 cos(x7) -sin(x7);
0 sin(x7)/cos(x8) cos(x7)/cos(x8)];

x7to9dot = H_phi*omega_b;

% Place in first order form

XDOT = [x1to3dot
x4to6dot
x7to9dot];

Interpreted MATLAB Function:

Pass the input values to a MATLAB function for evaluation. The function must return a single
value having the dimensions specified by 'Output dimensions' and 'Collapse 2-D results to 1-
D'.

Double-click the Non-linear 6DOF Aircraft Model shown in Fig. 8.1

90
 Fig. 8.2 Interpreted MATLAB Function Block Parameters

Fig. 8.3 Non-linear 6 DOF Model with limit for Input Values

Impulse
Input

Double-click on the Saturation Limit Double-click on the Control surface

Subsystem Model deflection Subsystem Model

91
Impulse Input Model:

Double-clicking the Step Model(s):

Example 8.2: Trim the aircraft model for steady-level flight condition.
Open the Simulink® model.

Fig. 8.6 Aircraft Model for Linearization (Trim Point Analysis)

92
Double-clicking the Subsystem block in the above model, we see

Fig. 8.7 Model Inside Subsystem Shown in Fig. 8.6

To open Model Linearizer for the model, in the Simulink model window, in the Apps gallery,
click Model Linearizer.

Fig. 8.8 Model Linearizer

Clicking the Model Linearizer

• To use the analysis points you defined in the Simulink model as linearization I/Os, on
the Linear Analysis tab, in the Analysis I/Os drop-down list, leave Model I/Os selected.

• To specify linearization input and output points, open the Linearization tab. To do so,
in the Apps gallery, click Linearization Manager.

• To specify an analysis point for a signal, click the signal in the model. Then, on the
Linearization tab, in the Insert Analysis Points gallery, select the type of analysis point.

▪ Configure the output signal of the Controller block as an Input Perturbation.

▪ Configure the output signal of the Magnetic Ball Plant block as an Open-loop
Output.

93
 Fig. 8.9 Linear Analysis Window

LINEAR ANALYSIS > Operating Point: > Trim Model

Fig. 8.10 Linear Analysis Window

➢ On selecting the Trim Model, a window pop-up as shown in Figure below to modify the
specifications.

➢ Click the States > Checkmark all the State Variables (State – 1, 2, 3, …).

➢ Then click the Outputs > Checkmark all the outputs and change the value of Output –
1 to 85.

➢ After the settings are changed > Click Start Trimming.

94
After clicking Start Trimming

95
Check the optimized results in the Linear Analysis Workspace

96
 Tutorial 9: Implementation of Kalman Filter in Simulink
Utilize the Kalman filter to enhance the estimation of the angular position theta in an inverted
pendulum system, considering the presence of noise in the measured values obtained from a
sensor, such as a rotary potentiometer.

Fig. 9.1 Inverted Pendulum

Input parameters:
Parameters Symbol Value

Mass of the cart (kg) M 0.5

Mass of the pendulum (kg) M 0.2
Friction of the cart (N/m/s) B 0.1
Length to pendulum center of mass (m) L 0.3
Inertia of the pendulum (kg*m^2) I 0.006
Force applied to the cart F ---
Cart position coordinate (m) X ---
Pendulum angle from vertical (radians) theta ---

Below are the free-body diagrams of the two elements of the inverted pendulum system.

Equations of Motion:

M x+ b x+ N = F

where,

97
 2
N = m x + ml θ cos θ − ml θ sin θ

Sum the moments about the centroid of the pendulum to get the following equation.

I θ = (1 I )( −Plsin θ − Nl cos θ )

where,

 2 
P = m  l θ cos ( θ ) + l θ sin ( θ ) + g 
 

Fig. 9.2 Model for Tutorial 9

Initial condition
Theta (0) = pi/18

Fig. 9.3 Model for Inverted Pendulum

98
 Fig. 9.4 Block Parameters

MATLAB Code:

clc;
clear;
close all;

%% Process Noise Covariance

Q = 0.0001;

%% Measurement Noise Covariance

R = 0.0005;

%% Sample Time

Ts = 0.01;

m = 0.5;
M = 0.2;
b = 0.1;
l = 0.3;
I = 0.006;
g = 9.81;

a11 = 0;
a12 = 1;
a13 = 0;
a14 = 0;
a21 = 0;
a22 = (-(I+m*l^2)*b)/(I*(M+m)+(M*m*l^2));
a23 = (m^2*g*l^2)/(I*(m+M)+M*m*l^2);
a24 = 0;
a31 = 0;
a32 = 0;
a33 = 0;

99
a34 = 1;
a41 = 0;
a42 = (-m*l*b)/(I*(M+m)+M*m*l^2);
a43 = ((m*g*l)*(M+m))/(I*(M+m)+M*m*l^2);
a44 = 0;

b21 = 0;
b22 = (I+m*l^2)/(I*(M+m)+M*m*l^2);
b23 = 0;
b24 = (m*l)/(I*(M+m)+M*m*l^2);

A = [a11 a12 a13 a14;a21 a22 a23 a24;a31 a32 a33 a34;a41 a42
a43 a44];
B = [b21;b22;b23;b24];
C = [1 0 0 0;0 0 1 0];
D = [0;0];

%% Run the Simulink Model

sim("Pend_Model.slx")

%% Plot the results

time = tout;

actual_theta = simout.Data(:,1);
noised_theta = simout.Data(:,2);
filtered_theta = simout.Data(:,3);
plot(time,actual_theta,'-r','linewidth',1.5)
hold on
plot(time,noised_theta,'--g','linewidth',1.5)
plot(time,filtered_theta,'-.b','linewidth',1.5)
xlabel('Time (sec)')
ylabel('Theta (rad)')
legend('Actual Theta','Noised Theta','Filtered Theta');

grid on
ax = gca;
ax.GridColor = [0 0 0];
ax.GridLineStyle = ':';
ax.GridAlpha = 0.3;
ax.Layer = 'bottom';

%% To save the figure in png format

print('Kalman Filter','-dpng','-r600');

100
 Fig. 9.5 Block Parameters

Fig. 9.6 Output from the Model

Introduction to Kalman Filter:

The Kalman filter keeps track of the estimated state of the system and the variance or
uncertainty of the estimate. The estimate is updated using a state transition model and

101
measurements. denotes the estimate of the system's state at time step k before the k-th
measurement yk has been taken into account; is the corresponding uncertainty.

Fig. 9.7 This diagram explains the basic steps of Kalman filtering: prediction and
update. It also illustrates how the filter keeps track of not only the mean value of the
state but also the estimated variance.

102
 Tutorial 10: Modeling and Simulation of Aerospace Vehicles
Open the Aerospace Blockset Library. You can access this library through the Simulink Library
Browser or directly open the Aerospace Blockset window from the MATLAB command line:

>> Aerolib

Fig. 10.1 Aerospace Blockset

If you prefer to open the complete model of a simple actuator system shown below instead of
building it, enter the following command at the MATLAB command window:

>> aeroblktutorial

To open the cockpit instruments, at the MATLAB command window, enter

>> aeroblk_HL20_Gauges

Fig. 10.2 HL20 Gauges/Instrument Panel

To view the airspeed correction models, enter the following at the MATLAB command line:

103
>> aeroblk_indicated
>> aeroblk_calibrated

Fig. 10.3 Airspeed Correction Model (TAS to IAS)

Fig. 10.4 Airspeed Correction Model (IAS to TAS)

To open the Wright Flyer model by entering at the MATLAB command line, enter

>> aeroblk_wf_3dof

104
Type the following in the MATLAB® command line. The model opens.

>> aeroblk_HL20

Fig. 10.6 HL-20 Space Shuttle Model

To open the Light Weight Airplane Design by entering at the MATLAB command line, enter

>> asbSkyHogg

Fig. 10.7 Light Weight Airplane Design Model

To open the Gravity Models with Precessing Reference Frame by entering at the MATLAB
command line, enter

>> asbGravWPrec

105
This model shows how to implement various gravity models with precessing reference frames
using Aerospace Blockset™ blocks. The Aerospace Blockset blocks are shown in red.

Fig. 10.8 Gravity Models with Precessing Reference Frame

To open the model used Airframe Trim and Linearize with Simulink Control Design by entering
at the MATLAB command line, enter

>> aeroblk_guidance_airframe

Fig. 10.9 Model Used in Airframe Trim and Linearization Model

106
Example 10.1: Thrust Vector Controlled Rocket Model
Equations of Motion:
Fx = T sin θ; Fz = T cos θ; M = Marm × Fx
where Fx = Force perpendicular to the rocket; Fz = Force along the longitudinal axis; M =
Moment about Center of Mass; Marm = Moment arm; Θ = Gimbal angle
Given Input:
Moment arm = 0.28 m; Mass = 0.543 kg; Mass moment of inertia = 0.048 kgm2

Block(s) Required:

• Step Input
• Slider with constant
• Degree to Radian
• Sin and cos block
• 3DOF (Body Axes)
• Gain
• XY Plot
• Scope(s)
• Compare to constant
• Stop

Fig. 10.10 Thrust Vector Controlled Rocket Model

Fig. 10.11 Input of the Rocket Model

107
 Fig. 10.12 Output of the Rocket Model
MATLAB Code:
clc;
clear;
close all;
TF = 11;
results = sim('Rocket_Model.slx');
plot(results.simout,'linewidth',2)
xlabel('Time (sec)')
ylabel('Height (m)')
title('Thrust Vector Controlled Rocket Model')
grid on
ax = gca;
ax.GridColor = [0 0 0];
ax.GridLineStyle = ':';
ax.GridAlpha = 0.3;
ax.Layer = 'bottom';

print('Rocket_Model','-dpng','-r600');

Fig. 10.13 Output for Example 10.1

108
Example 10.2: Implementation of MATLAB Animation Block (Create six-degrees-of-freedom
multibody custom geometry block.)
Aerospace Blockset > Animation > MATLAB-Based Animation > MATLAB Animation

Fig. 10.14 MATLAB Animation Model

MATLAB Code (Execute this code before simulate the Simulink model)

clc;
clear;
close all;

load simdata

Time = simdata(:,1);
X = simdata(:,[1,2]);
Y = simdata(:,[1,3]);
Z = simdata(:,[1,4]);
Roll = simdata(:,[1,5]);
Pitch = simdata(:,[1,6]);
Yaw = simdata(:,[1,7]);

Double-click the MATLAB Animation Block

Block Parameter

109
Block Parameter for MATLAB Animation Block

Fig. 10.15 Block Parameters

AC3D Models:

bluewedge.ac

delta2.ac

pa24–250_blue.ac

pa24–250_orange.ac

redwedge.ac

testrocket.ac

110
Example 10.3: Implementation of Simulink Animation Block

This part of model is taken from aeroblk_HL20

Fig. 10.16 Animation Model

Double-click on Visualization Subsystem to get the above part of the model.

Fig. 10.17 HL-20 Airframe Model

111
Example 10.4: Implementation of Simulink Animation Block

Fig. 10.18 Non-linear 6 DOF Aircraft Model with Navigation Equation

MATLAB Code for the model shown in Fig. 10.18

function V_NED = fcn(x1,x2,x3,x7,x8,x9)

C1v = [cos(x9) sin(x9) 0;-sin(x9) cos(x9) 0;0 0 1];

C21 = [cos(x8) 0 -sin(x8);0 1 0;sin(x8) 0 cos(x8)];
Cb2 = [1 0 0;0 cos(x7) sin(x7);0 -sin(x7) cos(x7)];

Cbv = Cb2*C21*C1v;
Cvb = Cbv';
V_b = [x1;x2;x3];

V_NED = Cvb*V_b;

Fig. 10.19 Model Inside Navigation Equation in Figure 10.18

112
 Fig. 10.20 Model Inside Visualization Subsystem

In-built Examples Using Aerospace Blockset(s)

Multiple Aircraft with Collaborative Control

>> asbswarm

Quaternion Estimate from Measured Rates

>> asbQuatEML

Modeling a Six Degree of Freedom Motion Platform

>> asb6Dof

Gravity Models with Precessing Reference Frame

>> asbGravWPrec

Airframe Trim and Linearize with Simulink Control Design

>> aeroblk_guidance_airframe

Quadcopter Project

>> asbQuadcopterStart

Electrical Component Analysis for Hybrid and Electric Aircraft

>> asbHybridAircraftStart

Constellation Modeling with the Orbit Propagator Block

>> OrbitPropagatorBlockExampleModel

Mission Analysis with the Orbit Propagator Block

113
>> OrbitPropagatorBlockExampleModel

Getting Started with the Spacecraft Dynamics Block

>> SpacecraftDynamicsBlockExampleModel

114