12

MATLAB Applications No. (1)

Modeling and Simulation

 Objectives:
The objectives of this application are to:
1. Convert the functions from time domain to Laplace domain (s).
2. Obtain the time domain function from the complex variable domain (s).
3. Solve the linear ordinary differential equations, and obtain the output
response.
4. Derive the mathematical model of mechanical, electrical, and fluid control
systems.
5. Simulate the derived mathematical model using MATLAB Simulink.
 Introduction:
To understand and control complex systems, one must obtain quantitative
mathematical models of these systems. It is necessary therefore to analyze the
relationships between the system variables and to obtain a mathematical
model. Because the systems under consideration are dynamic in nature, the
descriptive equations are usually differential equations. Furthermore. If these
equations can be linearized, then the Laplace transform can be used to simplify
the method of solution.

 Laplace Transform:

The Laplace transform method substitutes relatively easily solved algebraic

equations for the more difficult differential equations. This method can only be
applied to linear ordinary differential equations.
Laplace transforming the linear ordinary differential equations in term of the
independent variable (t) converts them into algebraic equations in the Laplace

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 13

transform variable (s). This provides a very convenient representation of

system dynamics.
The operation of Laplace transformation will be indicated by the notation:
∞
ℒ[𝑓(𝑡)] = ∫ 𝑓(𝑡) 𝑒 −𝑠𝑡 𝑑𝑡 = 𝐹(𝑠).
0

Where:
ℒ: Laplace transform operator.
s: Laplace transform variable (this variable is a complex number).
Example 2.1: Use MATLAB to find the Laplace transform of the following
functions:

1. f(t) = 𝑢(𝑡). 𝟐. f(t) = 𝑡 𝑢(𝑡). 𝟑. f(t) = 𝑒 −𝑎𝑡 .

4. f(t) = sin(𝜔𝑡). 5. f(t) = cos(𝜔𝑡).

Solution:
>> % 1. To determine the Laplace transform for unit step function %
>> syms t
>> ft=heaviside(t);
>> Fs=laplace(ft);
>> pretty(Fs)
1
------
s
>> % 2. To determine the Laplace transform for unit ramp function %
>> syms t
>> ft=t;
>> Fs=laplace(ft);
>> pretty(Fs)
1
------
s2

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 14

>> % 3. To determine the Laplace transform for exponential function %

>> syms t a
>> ft=exp(-a*t);
>> Fs=laplace(ft);
>> pretty(Fs)
1
-------
s+a

>> % 4. To determine the Laplace transform for sine function %

>> syms t w
>> ft=sin(w*t);
>> Fs=laplace(ft);
>> pretty(Fs)
w
----------
s2 + w2

>> % 5. To determine the Laplace transform for cosine function %

>> syms t w
>> ft=cos(w*t);
>> Fs=laplace(ft);
>> pretty(Fs)
s
----------
s2 + w2

Example 2.2: Use MATLAB to find the Laplace transform of f(t) are defined by
the following function:

1. f(t) = 5𝑒 −2𝑡 cos(5𝑡).

2. f(t) = (𝑡 + 1)2 + 7 𝑢(𝑡).

3. f(t) = 𝑒 𝑡 + 2 𝑒 −2𝑡 + 𝑡 𝑒 −3𝑡 .

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 15

Solution:
>> % 1. To find the Laplace transform of f(t)= 5 exp(-2t) cos(5t) %
>> syms t
>> ft=5*exp(-2*t)*cos(5*t);
>> Fs=laplace(ft);
>> pretty(Fs)
5 (s + 2)
------------------
(s + 2)2 + 25

>> % 2. To find the Laplace transform of f(t)=(t+1)^2+7u(t) %

>> syms t
>> ft=(t+1)^2+7;
>> Fs=laplace(ft);
>> Fs=simplify(Fs);
>> pretty(Fs)

2 (4 s2 + s + 1)
----------------------
s3
>> % 3. To find the Laplace transform of f(t)=exp(t)+2 exp(-2t) + t exp(-3t)%
>> syms t
>> ft=exp(t)+2*exp(-2*t)+t*exp(-3*t);
>> Fs=laplace(ft);
>> pretty(Fs)
1 2 1
----- + ----- + ----------
s - 1 s + 2 (s + 3)2

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 16

 Inverse Laplace Transform:

After transforming equations into the Laplace domain and solving for output
variables as functions of (s), we sometimes want transform back into the time
domain. This equation is called inversion of Laplace transforms. The inverse
Laplace transform is written as:

𝜎+𝑗∞
1
𝑓(𝑡) = ∫ 𝐹(𝑠) 𝑒 𝑠𝑡 𝑑𝑠.
2𝜋𝑗
𝜎−𝑗∞

The inverse Laplace transformation is usually obtained by using partial fraction

expansion. The function to be inverted 𝐹(𝑠) is merely rearranged into a series
of simple function:

𝐹(𝑠) = 𝐹1 (𝑠) + 𝐹2 (𝑠) + ⋯ + 𝐹𝑁 (𝑠)

Then each term is inverted (usually by inspection because they are simple). The
total time dependent function in the sum of the individual time dependent
functions:

𝑓(𝑡) = ℒ −1 [𝐹1 (𝑠)] + ℒ −1 [𝐹2 (𝑠)] + ⋯ + ℒ −1 [𝐹𝑁 (𝑠)]

𝑓(𝑡) = 𝑓1 (𝑡) + 𝑓2 (𝑡) + ⋯ + 𝑓𝑁 (𝑡)

Example 2.3: Find the inverse Laplace transform of the following functions.
First, perform the partial fraction expansion on Y(s).
4
𝟏. Y(𝑠) =
𝑠(𝑠 + 2)(𝑠 + 3)
(𝑠 + 5)
𝟐. Y(𝑠) =
(𝑠 + 4)2 (𝑠 + 7)

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 17

Solution:

>> % To perform the partial fraction expansion for: y(s)=4/(s(s+2)(s+3))%

>> num=[4];
>> den1=[1 0];den2=[1 2];den3=[1 3];
>> den=conv(den1,conv(den2,den3));
>> [r,p,k]=residue(num,den)
r=
1.3333
-2.0000
0.6667
p=
-3.0000
-2.0000
0
k=
[]
>> % To find the inverse Laplace transform …%

>> syms s
>> Ys=1.33/(s+3)-2/(s+2)+0.6667/s;
>> yt=ilaplace(Ys)
yt =
1.33 * exp(-3t) - 2 * exp(-2t) +0.6667

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 18

>> % To perform the partial fraction expansion for: y(s)=(s+5)/((s+4)2(s+7))%

>> num=[1 5];
>> den1=[1 4];den2=[1 4];den3=[1 7];
>> den=conv(den1,conv(den2,den3));
>> [r,p,k]=residue(num,den)
r=
-0.2222
0.2222
0.3333
p=
-7.0000
-4.0000
-4.0000
k=
[]

>> % To find the inverse Laplace transform …%

>> syms s
>> Ys=-0.22/(s+7)+0.22/(s+4)+0.33/(s+4)^2;
>> yt=ilaplace(Ys)
yt =
0.22*exp(-4*t) – 0.22*exp(-7*t) + 0.33*t *exp(-4*t)

Example 2.4: Use MATLAB to find the inverse Laplace transform directly of the
following functions:

1
𝟏. Y(𝑠) =
(𝑠 + 1)(𝑠 2 + 4)
𝑠 2 + 6𝑠 + 8
𝟐. Y(𝑠) = 2
𝑠 + 8𝑠 + 15

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 19

Solution:
>> % To find inverse Laplace transform for: y(s)=1/((s+1)(s^2+4)) directly %
>> syms s
>> Ys=1/((s+1)*(s^2+4));
>> yt=ilaplace(Ys)
yt =
1/5*exp(-t) - cos(2*t)/5 + sin(2*t)/10
>> % To find inverse Laplace transform for…%
>> % y(s)=(s^2+6s+8)/(s^2+8s+15) directly %
>> syms s
>> Ys=(s^2+6*s+8)/(s^2+8*s+15);
>> yt=ilaplace(Ys)
yt =
dirac(t) - 3/2*exp(-5*t) - 1/2*exp(-3*t)

Example 2.5: An RC Lag circuit is shown in Figure 2.1. The output of the circuit
is described by the following expression in the Laplace domain:
10
𝑉𝑜 (𝑠) =
𝑠(𝑠 + 10)

Figure 2.1: An RC Lag circuit.

1. Use the inverse Laplace transform to find 𝑉𝑜 (𝑡).
2. Plot the response of the output.
3. Repeat parts (1) and (2) using MATLAB Simulink.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 20

Solution:
>> % 1. To find inverse Laplace transform for: Vo(s)=10/(s(s+10))%
>> syms s;
>> Vs=10/(s*(s+10));
>> vt=ilaplace(Vs)
vt =
1 - exp(-10*t)
>> % 2. To plot output response of the RC Lag network: v(t)= 1-exp(-10*t)%
>> t=0:0.01:1;
>> vt=1-exp(-10*t);
>> plot(t,vt)
>> grid
>> title('Output Response of: v(t)= 1-exp(-10*t)')
>> xlabel('Time (sec)')
>> ylabel('Output Vo(Voltage)')

Output Response of: v(t)= 1-exp(-10*t)

1

0.9

0.8

0.7
Output Vo(Voltage)

0.6

0.5

0.4

0.3

0.2

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time (sec)

Figure 2.2: Output response of RC Lag network for unit step input.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 21

3. From Figure 2.1. The equations for RC circuit are:

𝑣𝑖𝑛 (𝑡) − 𝑣𝑜 (𝑡)
𝑖(𝑡) = (2.1)
𝑅
1
𝑣𝑜 (𝑡) = ∫ 𝑖(𝑡) 𝑑𝑡 (2.2)
𝐶
From Equations (2.1) and (2.2), we can obtain the elements of the block
diagram, as shown in Figure 2.3(a).

By connecting signals properly, we can construct a block diagram, as shown in

Figure 2.3(b).

Figure 2.3: (a) Elements of the block diagram of the RC network shown in
Figure 2.1, (b) block diagram of the RC network.
The Simulink diagram of the coupled tank system as shown in Figure 2.4:

Figure 2.4: Simulink diagram of the RC circuit.

The output voltage response that has been obtained from the simulation is the
same as the output response of RC Lag network for unit step input shown in
Figure 2.2.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 22

Example 2.6: Consider the mass m mounted on a massless cart, as shown in

Figure 2.5. The external force f(t) is the input to the system (f(t) = 𝑒 −3𝑡 ), and
the displacement 𝑥(𝑡) of the mass is the output. The force equation of motion
in the 𝑥 direction:

𝑓(𝑡) = 𝑚𝑥̈ + 𝑏𝑥̇ + 𝑘𝑥

If:m = 5 Kg, 𝑏 = 15(N. sec/m) , 𝑘 = 10 N/m, and the initial conditions (𝑥(0) =
0.5 m, 𝑥̇ (0) = 3 m/sec).

Using MATLAB:
1. Solve the differential equation above to determine 𝑥(𝑡).
2. Plot the output response of the system.

Figure 2.5: Mass on cart.

Solution:

>> % 1. To find inverse Laplace transform to solve DE 𝐟(𝐭) = 𝐦𝐱̈ + 𝐛𝐱̇ + 𝐤𝐱 %

>>% When f(t)=exp(-3*t), and the initial conditions: x(0)=0.5m, 𝐱̇ (𝟎) =3 m/sec.%

>> xt=dsolve('5*D2x+15*Dx+10*x=exp(-3*t)','x(0)=0.5,Dx(0)=3')
xt =
4.1*exp(-t) – 3.7*exp(-2*t) + 0.1*exp(-3*t)

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 23

>> % 2. To plot the output response of …%

>> % x(t)=4.1*exp(-t)-3.7*exp(-2*t)+0.1*exp(-3*t) %
>> t=0:0.1:10;
>> xt=4.1*exp(-t)-3.7*exp(-2*t)+0.1*exp(-3*t);
>> plot(t,xt)
>> grid
>> title('The Output Response for f(t)=exp(-3t)')
>> xlabel('Time (sec)')
>> ylabel('Diplacement (Meter)')

The Output Response for f(t)=exp(-3t)

1.4

1.2

1
Diplacement (Meter)

0.8

0.6

0.4

0.2

0
0 1 2 3 4 5 6 7 8 9 10
Time (sec)
Figure 2.6: Output response of spring mass damper system for exponential input.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 24

Example 2.7: Consider the mechanical system shown in Figure 2.7. The external
force f(t) is the input to the system, and the displacement y(t) of the mass is the
output. This system is a single-input, single-output system:
𝑚𝑦̈ + 𝑏𝑦̇ + 𝑘𝑦 = 𝑓(𝑡)

Figure 2.7: Mechanical system.

Simulate the system and obtain the output response using MATLAB Simulink.

Solution:
The simulation of the system using MATLAB Simulink as shown in Figure 2.8:

Figure 2.8: Simulation of the mechanical system of Example 2.7 using MATLAB Simulink.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 25

The step response of the mechanical system in Example 2.7 as shown in Figure
2.9:

Output Response of Mechanical System

1.5

1
Diplacement(Meter)

0.5

0
0 1 2 3 4 5 6 7 8 9 10
Time(sec)

Figure 2.9: Step response of the mechanical system for Example 2.7.

Example 2.8: A coupled spring – mass – damper system is shown in Figure

2.10. This system has one input which is 𝑢(𝑡)and two outputs which are the
displacement 𝑥1 and 𝑥2 of the first and second cart respectively.

Figure 2.10: Two cart mechanical system.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 26

From the Figure 2.10, the systems equations at zero initial conditions are:
𝑚1 𝑥1̈ + 𝑏(𝑥1̇ − 𝑥2̇ ) + 𝑘1 𝑥1 + 𝑘2 (𝑥1 − 𝑥2 ) = 𝑢(𝑡)
𝑚2 𝑥2̈ + 𝑏(𝑥2̇ − 𝑥1̇ ) + 𝑘3 𝑥2 + 𝑘2 (𝑥2 − 𝑥1 ) = 0
Use MATLAB Simulink to:
1. Simulate the two cart mechanical system.
2. Obtain the response of 𝑥1 with respect to 𝑢(𝑡).
Solution:
The simulation of the two cart mechanical system using MATLAB Simulink as
shown in Figure 2.11:

Figure 2.11: Simulink diagram of the two cart mechanical system.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 27

The sinusoidal response of 𝑥1 with respect to 𝑢(𝑡) of the two cart mechanical
system as shown in Figure 2.12:
Output Response of Mechanical System
0.015

0.01

0.005
Diplacement x1(Meter)

-0.005

-0.01

-0.015
0 5 10 15 20 25 30 35 40 45 50
Time(sec)

Figure 2.12: Sinusoidal response of the two cart mechanical system.

Example 2.9: Consider the system shown in Figure 2.13. In this system, the two
tanks interact. Thus, the transfer function of the system is not the product of
two first-order transfer functions.

Figure 2.13: Coupled tank control system.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 28

Using the symbols as defined in Figure 2.13, we can derive the following
equations for this system:
ℎ1 − ℎ2
= 𝑞1 (2.3)
𝑅1
𝑑ℎ1
𝐶1 = 𝑞 − 𝑞1 (2.4)
𝑑𝑡
ℎ2
= 𝑞2 (2.5)
𝑅2
𝑑ℎ2
𝐶2 = 𝑞1 − 𝑞2 (2.6)
𝑑𝑡
Simulate the system using MATLAB Simulink to obtain the step response of
the outlet flow rate.
Solution:
From Equations (2.3) through (2.6), we can obtain the elements of the block
diagram, as shown in Figure 2.14(a).

By connecting signals properly, we can construct a block diagram, as shown in

Figure 2.14(b).

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 29

Figure 2.14: (a) Elements of the block diagram of the system shown in
Figure 2.10; (b) block diagram of the system.

The Simulink diagram of the coupled tank system as shown in Figure 2.15:

Figure 2.15: Simulink diagram of the coupled tank system.

The step response of the coupled tank control system as shown in Figure 2.16:
Output Response of Liquid Flow System

1
Outlet flow rate(m3/sec)

0.8

0.6

0.4

0.2

0
0 5 10 15 20 25 30 35 40 45 50
Time(sec)

Figure 2.16: Step response of the coupled tank system.

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 30

Homework 2.1: Find the Laplace transform of the following time functions
using MATLAB:
1. f(t) = 𝜋 + 7𝑡 + 𝑡 2 + 3 𝛿(𝑡) (Where 𝛿(𝑡): Impulse function).

2. f(t) = sinh 2𝑡 + 2 cos 2𝑡 + 𝑒 −2𝑡 sin 2𝑡

3. f(t) = sin2 𝑡 + 3 cos 2 𝑡

---------------------------------------------------------------------------------------------------
Homework 2.2: Use MATLAB to determine the inverse Laplace transform of the
following transfer functions:
𝑒 −2𝑠
𝟏. Y(𝑠) = (Directly)
𝑠(𝑠 + 4)(𝑠 + 6)
2𝑠 4 + 9𝑠 3 + 15𝑠 2 + 𝑠 + 2
𝟐. Y(𝑠) = (Using Residue Command)
𝑠 2 (𝑠 + 2)(𝑠 + 1)2
𝑠 3 + 7𝑠 2 + 9𝑠 + 5
𝟑. Y(𝑠) = (Using Residue Command)
(𝑠 + 3)(𝑠 + 4)
---------------------------------------------------------------------------------------------------
Homework 2.3: A ship has a mass 𝑚 and resistance 𝑐 times the forward
velocity 𝑢(𝑡). If the thrust from the propeller is 𝑘 times its angular velocity
𝜔(𝑡), where:
𝑑𝑢𝑓
𝑚 + 𝑐𝑢𝑓 (𝑡) = 𝑘𝜔(𝑡)
𝑑𝑡
If: 𝑚 = 180 × 105 kg, 𝑐 = 150 × 103 𝑁. 𝑠/𝑚, 𝑘 = 96 × 103 𝑁. 𝑠/rad , and
𝜔(𝑡) = 12.5 𝑢(𝑡).
Find the following using MATLAB:
1. Solve the differential equation above to determine 𝑢𝑓 (𝑡).
2. Plot the output response of a ship system.
3. From part (2), what is the forward velocity after one minute and after ten
minutes?

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 31

Homework 2.4: The rotational velocity ω of a satellite shown in Figure 2.17 is

adjusted by changing the length L of a beam (𝐿(𝑠) = 1/4𝑠 ). The rotational
velocity ω is described by the following expression in the Laplace domain:

2.5(𝑠 + 2)
𝜔(𝑠) = 𝐿(𝑠)
(𝑠 + 5)(𝑠 + 1)2

Figure 2.17: Satellite with adjustable rotational velocity.

1. Perform the partial fraction expansion on 𝜔(𝑠).
2. Use the inverse Laplace transform to find 𝜔(𝑡).
3. Plot the output response rotational velocity of the satellite.
-----------------------------------------------------------------------------------------------------
Homework 2.5: A load added to a truck results in a force F on the support spring,
and the tire flexes as shown in Figure 2.18(a). The model for the tire movement
is shown in Figure 2.18(b).
Use MATLAB Simulink to:
1. Simulate the truck support mechanical system.
2. Obtain the unit step response of 𝑥1 and 𝑥2 with respect to f(t).

Khaled Mustafa Mahmoud Session: Spring 2017/2018

 32

Figure 2.18: Truck support model.

-------------------------------------------------------------------------------------------------
Homework 2.6: For a two series tank control system shown in Figure 2.19.
The inlet flow rate 𝑞𝑖 is manipulated using the control valve to regulate the
outlet flow rate 𝑞0 .

Figure 2.19: Series tank control system.

Using MATLAB Simulink:
1. Simulate the system model.
2. Obtain the step response of the outlet flow rate.

Khaled Mustafa Mahmoud Session: Spring 2017/2018