EEE116 Experimental, Computer Skills and Sustainability: MATLAB Group Assignment
Due on Monday, 13th May, 2019, 18:00
1 / 16
Assignment Regulations
⚫ This is a group assignment. Each group MUST submit one soft copy of the group report via the ICE before the due date.
⚫ A coversheet can be created in your own way, but the following information MUST be included: your group number, student ID number, full name and email address of each group member.
⚫ All the formula, derivations, completed MATLAB scripts and functions with original highlighted text format in MATLAB editor, computational results in the command window, and plotted figures, should be part of the answers.
⚫ There is no hard requirement on how the report must be organized. You can organize your report question by question (i.e. give one section for each question). Then, for each section, you can organize it in your own way. However, the contents and information required by each individual question MUST be provided.
⚫ You may refer to textbooks and lecture notes to discover approaches to problems, however, the assignment should be your own group work.
⚫ Where you do make use of other reference, please cite them in your work. Students are reminded to refer and adhere to plagiarism policy and regulations set by XJTLU. References, in IEEE style can be attached as an appendix.
⚫ A contribution form MUST be provided which can be given as the final page of your report. You can follow the example contribution form provided in the final page of this document.
⚫ Assignments may be accepted up to 5 working days after the deadline has passed; a late penalty of 5% will be applied for each working day late without an extension being granted. Submissions over 5 working days late will not be marked. Emailed submissions will NOT be accepted without exceptional circumstances.
2 / 16
MATLAB Group Assignment for EEE116
Part I: developing basic programming skills by solving simple engineering problems
Question I-A, Numerical method for solution of ordinary differential equations. (25 Marks)
The classical method of solving ordinary differential equations (ODE) is the 4th order Rugge-Kutta method. The algorithm is given below:
Assume we have a 1st order ODE:
𝑑𝑦 = 𝑓(𝑡, 𝑦) 𝑑𝑡
For the given time period 𝑡 ∈ [𝑡0, 𝑡𝑓], we need to find out the time variation of y, given
the initial condition 𝑦(𝑡0) = 𝑦0. This kind of problem is referred to as initial value problems.
For solution, let us have divide the time period into N intervals (i.e. N+1 time instants), then the value of y at (n+1)th time instant is given by the following relations with n from 0 to N with h as the time interval size:
(I-A.1)
where
𝑦 =𝑦 +h(𝑘 +2𝑘 +2𝑘 +𝑘) 𝑛+1 𝑛 6 1 2 3 4
𝑘 =𝑓(𝑡,𝑦) 1𝑛𝑛
𝑘 =𝑓(𝑡 +h,𝑦 +h𝑘) 2𝑛2𝑛21
𝑘 =𝑓(𝑡 +h,𝑦 +h𝑘) 3𝑛2𝑛22
𝑘 =𝑓(𝑡 +h,𝑦 +h𝑘) 4𝑛𝑛3
(I-A.2)
(I-A.3)
Question (a) (15 Marks): Now, solve an RC circuit problem: the equation for the voltage y across the capacitor is
𝑅𝐶 𝑑𝑦 + 𝑦 = 𝑣(𝑡) 𝑑𝑡
(I-A.4) where 𝑣(𝑡) is the applied voltage to the circuit. Suppose that RC = 0.2 s and that the capacitor voltage, y, is initially 2 V (i.e. 𝑦(0) = 2 𝑉). Suppose also that the applied
3 / 16
voltage is 𝑣(𝑡) = 10[2 − 𝑒−𝑡 sin(5𝜋𝑡)] V. Plot the voltage y(t) for 0 ≤ 𝑡 ≤ 5 s. Requirements:
(1) (7 Marks) Write a function file to give a general ODE solver using Runge-Kutta method, e.g.
function [ t,y ] = marunge4s( dyfun,tspan,y0,h )
where “marunge4s” is the function name; “dyfun” is the name of the ODE to be solved which is defined by using a function; “tspan” is the time period to be simulated; y0 is the initial condition of the solved-for variable, “h” is the time interval size.
(2) (5 Marks) Do the above question. i.e. solve Equation I-A.4, by calling the function created. Then, plot the voltage y(t) for 0 ≤ 𝑡 ≤ 5 s.
Hints: For this question, you need to create your own function file using Runge-Kutta method (This is done for Question (1)), which is similar to ode45 (a MATLAB standard function of solving ODEs), to solve the ODE defined by another function (refer to hints of Question (3) for how to do this). Then, the calling pattern of your function (to solve ODEs) should be, e.g.
[t,y1]=marunge4s(‘RC_circuit’,[0 0.5],2,0.001);
where “RC_circuit” is the name of the ODE you have defined using a function file; 0.001 is h, the time interval size. You need to use the function ‘feval’ in your function “marunge4s” if the pattern is shown above. If you do not wish to use ‘feval’, then the calling pattern should be:
[t,y1]=marunge4s(@RC_circuit,[0 0.5],2,0.001);
(3) (3 Marks) Check your solution by using a MATLAB standard function for solving ODEs, i.e. ode45: you need to plot results of voltage y(t) for 0 ≤ 𝑡 ≤ 5 s, obtained by using both your user defined function and ode 45, then, compare the results and comment on the features (e.g. whether they are the same or not?)
Hints: An example of using ode45: we need to first define the ODE to be solved which can be done by using a function, e.g.
The function file defining the ODE to be solved:
function [ ydot ] = RC_circuit( t,y )
%define the ODE by using a function to be called by
%MATLAB function ode45 or
%your own function solver using Runge-Kutta method
Vsource=100*sin(2*pi*50*t);
ydot=-5*y+5*Vsource;
end
4 / 16
Then, use ode45 function to solve the equation you just defined using a function, as shown below:
The main programme written in a script file:
clear all
close all
clc
%
%’RC_circuit’ is the ODE name defined by a function
%[0 0.5] is the time period for the simulation
%2 is the initial conditions, v_c(t=0)=2 V
[t,y]=ode45(‘RC_circuit’,[0 0.5],2);
plot(t,y(:,1))
xlabel(‘time (s)’)
ylabel(‘capacitor voltage (V)’)
Question (b) (10 Marks): Now, we consider the RLC circuit with R, L and C connected in series, in which the capacitor voltage vc is governed by a 2nd order ODE:
𝐿𝐶 𝑑2𝑣𝑐 + 𝑅𝐶 𝑑𝑣𝑐 + 𝑣𝑐 = 𝑣(𝑡) 𝑑𝑡2 𝑑𝑡
(I-A.5) where v(t) is the applied power supply voltage. Suppose we have R = 100 ohms, L =
0.5 Henries and C = 2e-6 Farads. Solve the equations with the following conditions: (1) (5 Marks) Model input: 𝑣(𝑡) = 0 𝑉; Initial conditions: 𝑣𝑐(0) = 100 𝑉, 𝑑𝑣𝑐 (0) =
𝑑𝑡
0, then plot the results, i.e. the voltage 𝑣𝑐(𝑡) for 0 ≤ 𝑡 ≤ 0.5 s.
(2) (5 Marks) Model input: 𝑣(𝑡) = 100 sin(2𝜋𝑓𝑡) V , with f = 50 Hz; Initial
conditions: 𝑣𝑐(0) = 100 𝑉 𝑎𝑛𝑑 0 𝑉, 𝑑𝑣𝑐 (0) = 0, then plot the results for both two 𝑑𝑡
initial conditions, i.e. the voltage 𝑣𝑐(𝑡) for 0 ≤ 𝑡 ≤ 0.5 s.
Hints: The 2nd order ODE can be written into a set of 2 1st oder ODEs (see below), and
then the question becomes to solve 2 1st order ODEs using Kugge-Kutta method. Assume two vector spaces, X and Y, with the size of 2 × 1 each, then we let:
x(1) = 𝑣𝑐 and x(2) = 𝑑𝑣𝑐 𝑑𝑡
xdot(1) = x(2) = 𝑑𝑣𝑐 and xdot(2) = 𝑑2𝑣𝑐
𝑑𝑡
Thus, the 2nd order ODE becomes two 1st order ODEs
𝑑𝑡2
5 / 16
𝑥𝑑𝑜𝑡(1) 𝑑𝑿
𝑓 (𝑡, 𝑥(1), 𝑥(2)) =[ 𝑑𝑡 ]=𝐹(𝑡,𝑿)=[1 ]
𝑑𝑥(1)
𝑿𝒅𝒐𝒕=[ ]= 𝑥𝑑𝑜𝑡(2)
𝑑𝑥(2) 𝑓 (𝑡, 𝑥(1), 𝑥(2)) 2
𝑑𝑡
with the initial conditions 𝑥(1)|𝑡=0 = 𝑥0(1) and 𝑥(2)|𝑡=0 = 𝑥0(2). Equation (I-A.5)
can now be written as
𝑑𝑡
𝑥𝑑𝑜𝑡(1) = 𝑑𝑥(1) = x(2) 𝑑𝑡
𝑥𝑑𝑜𝑡(2)= 1 𝑣(𝑡)−𝑅𝑥(2)− 1 𝑥(1) 𝐿𝐶 𝐿 𝐿𝐶
(I-A.6)
(I-A.7)
with initial conditions 𝑥(1)|𝑡=0 = 𝑣𝑐|𝑡=0 and 𝑥(2)|𝑡=0 = 𝑑𝑣𝑐|
𝑑𝑡 𝑡=0
The problem then becomes to solve the vector equation:
𝑑𝑿 = 𝐹(𝑡,𝑋) 𝑜𝑟 𝑿𝒅𝒐𝒕 = 𝐹(𝑡,𝑋) 𝑑𝑡
(I-A.8) with 𝑿𝒅𝒐𝒕 = [𝑥𝑑𝑜𝑡(1)] and 𝑿 = [𝑥(1)]. For using MATLAB supplied function,
𝑥𝑑𝑜𝑡(2) 𝑥(2)
ode45, an example is given below:
The function file defining the ODE to be solved:
function [ xdot ] = CircuitRLC( t,x )
%define the ODE by using a function to be called by %MATLAB function ode45 or
%your own function solver using Runge-Kutta method Vsource=100*sin(2*pi*50*t);
C=2e-6;
R=100;
L=0.5;
xdot=zeros(2,1);
xdot(1)=x(2); xdot(2)=Vsource/(L*C)-(R/L)*x(2)-(1/(L*C))*x(1); end
The main programme written in a script file:
clear all
close all
clc
%
6 / 16
%[0 0.5] is the time period for the simulation
%[100 0] is the initial conditions, v_c(t=0)=100 V,
dv_c/dt(t=0)=0
[t,x]=ode45(‘CircuitRLC’,[0 0.3],[100 0]);
plot(t,x(:,1))
xlabel(‘time (s)’)
ylabel(‘capacitor voltage (V)’)
The user defined function file required in Question (a), “marunge4”, should also be able to solve this kind of 2nd order ODEs, not just one 1st order ODE. Think about how to do this.
Requirements: for Question I-A, in your report, the following information needs to be included:
a) Circuit diagram, the relevant equation derivations and all the known conditions stated in the question;
b) Flow chart of the programme design;
c) Coding with detailed notes on the key steps;
d) Results.
7 / 16
Question I-B, Kinematics analysis for a slider-crank linkage.
(25 Marks)
Consider a slider-crank linkage which is the fundamental mechanism used in car engines, a schematic diagram is shown in Figure I-B.1. Knowing the geometrical parameters of the mechanism and the rotational speed of the driver linkage, try to do the kinematics analysis for the mechanism to calculate the displacements, velocity and acceleration.
Figure I-B.1. Schematic of a slider-crank linkage mechanism. The general procedure of kinematic analysis is given below:
(a) Displacement analysis
The position vector equation is
⃗ ⃗ ⃗⃗⃗ 𝐿1 + 𝐿2 = 𝑆𝑐
Write it in the form of complex numbers, it becomes
𝐿1𝑒𝑖𝜃1 +𝐿2𝑒𝑖𝜃2 =𝑆𝑐
Separate the real and imaginary parts, we have
(I-B.1) (I-B.2)
(I-B.3)
(I-B.4)
Then we have
𝐿1 cos 𝜃1 + 𝐿2 cos 𝜃2 = 𝑆𝑐 𝐿1 sin 𝜃1 + 𝐿2 sin 𝜃2 = 0
𝜃2 = arcsin (−𝐿1 sin 𝜃1) 𝐿2
𝑆𝑐 = 𝐿1 cos 𝜃1 + 𝐿2 cos 𝜃2
8 / 16
(b) Velocity analysis
Find the first derivative of equation (B.2), we have
i𝐿1𝜔1𝐿1𝑒𝑖𝜃1 + i𝐿2𝜔2𝐿1𝑒𝑖𝜃2 = 𝑣𝑐
Separate the real and imaginary parts, we have
𝐿1𝜔1 cos 𝜃1 + 𝐿2 𝜔2cos 𝜃2 = 0 −𝐿1𝜔1sin𝜃1 −𝐿2𝜔2sin𝜃2 =𝑣𝑐
Write the above equation in the matrix form
[ 𝐿2 sin𝜃2 1][𝜔2] = 𝜔1 [−𝐿1 sin𝜃1] −𝐿2cos𝜃2 0 𝑣𝑐 𝐿1cos𝜃1
(c) Acceleration analysis
(I-B.5)
(I-B.6)
(I-B.7)
Get the second derivative for equation (B.2) then separate real and imaginary parts (equivalent to directly differentiate equation (B.6)), we have
𝐿 𝛼 cos𝜃 −𝐿 𝜔2sin𝜃 +𝐿 𝛼 cos𝜃 −𝐿 𝜔2sin𝜃 =0 111111222122
−𝐿 𝛼 sin𝜃 −𝐿 𝜔2cos𝜃 −𝐿 𝛼 sin𝜃 −𝐿 𝜔2cos𝜃 =𝑎 111111222222𝑐
Write the above equation in the matrix form
[ 𝐿2sin𝜃2 1][𝛼2]+[𝜔2𝐿2cos𝜃2 0][𝜔2] −𝐿2cos𝜃2 0 𝑎𝑐 𝜔2𝐿2sin𝜃2 0 𝑣𝑐
= 𝜔1 [−𝜔1𝐿1 cos 𝜃1] + 𝛼1 [−𝐿1 sin 𝜃1] −𝜔1𝐿1 sin 𝜃1 𝐿1 cos 𝜃1
(I-B.8)
(I-B.9)
For the mechanism of Figure I-B.1, suppose L1 = 100 mm, L2 = 300 mm, 𝜔1 = 10 𝑟𝑎𝑑/𝑠, the dimensions of the slider block is [x,y]=[80 mm,20 mm], try to do the following:
(1) (10 Marks) Calculate the displacement, velocity and acceleration for Slider C (sc, vc, ac) and the angular counterparts for Linkage 2 (𝜃2, 𝜔2, 𝛼2). You can create a function file (e.g. [outputs] = slider_crank(inputs)) for this by using Equations (I-B.4), (I-B,7) and (I-B.9).
(2) (5 Marks) Plot the results as a function of the driver rotational angle (𝜃1) which ranges from 0 to 720 deg. Call the function created in Question (1) to get the results then plot them out.
(3) (10 Marks) Create the movie showing the movement of this mechanism, when the driver rotational angle (𝜃1) varies from 1 to 360 deg. You may need to use the MATLAB functions “getframe” and “movie()”.
9 / 16
Requirements: for Question I-B, in your report, the following information needs to be included:
a) The schematic diagram of the mechanism, the relevant equation derivations and all the known conditions stated in the question;
b) Flow chart of the programme design;
c) Coding with detailed notes on the key steps;
d) Results.
10 / 16
Problem I-C, Solve nonlinear equations using Newton-Raphson method.
(10 Marks)
Write a programme to solve the circuit with a diode as shown in Figure I-C.1.
Figure I-C.1. A simple diode circuit, where VPS = 5 V and R = 2 kohms. The current-voltage characteristic of the diode is
(𝑉𝐷 ) 𝐼𝐷=𝐼𝑠[𝑒𝑉𝑇 −1]
where 𝐼 = 10−13 𝐴, and, 𝑉 = 0.026 𝑉. 𝑠𝑇
The governing equation of the circuit yields:
(I-C.1)
(𝑉𝐷 )
𝑉 =𝐼𝑅[𝑒𝑉𝑇 −1]+𝑉 𝑃𝑆𝑠 𝐷
(I-C.2) Equation (I-C.2) is nonlinear which is difficult to be solved analytically. Therefore, try to write a program using Newton-Raphson algorithm to solve it, then print the result of the solved value VD in the command window. Details of the Newton-Raphson method and the example has ben given in MATLAB lecture notes “Matlab Tutorial 2-2018-19”,
Slides 50-55.
Requirements: for Question I-C, in your report, the following information needs to be included:
a) Circuit diagram, the relevant equation derivations and all the known conditions stated in the question;
b) Flow chart of the programme design;
c) Coding with detailed notes on the key steps;
d) Results.
11 / 16
Part II: Use Simulink to solve simple engineering problems Question II-A, The motion of a two-mass suspension model.
(20 Marks)
Consider the two-mass suspension model shown in Figure II-A.1. The equations of motion are:
𝑚1𝑥̈1 +(𝑐1 +𝑐2)𝑥̇1 +(𝑘1 +𝑘2)𝑥1 −𝑐2𝑥̇2 −𝑘2𝑥2 =0 𝑚2𝑥̈2 + 𝑐2𝑥̇2 + 𝑘2𝑥2 − 𝑐2𝑥̇1 − 𝑘2𝑥1 = 𝑓(𝑡)
Supposethat 𝑚1 =𝑚2 =1, 𝑐1 =3, 𝑐2 =1, 𝑘1 =1 and 𝑘2 =4.
(II-A.1)
(1) (15 Marks) Develop a Simulink model of this system. In your report, you need to give information on:
a) Answer a question: whether to use a state-variable representation or a transfer-
function representation of the model, and why?
(3 Marks)
b) How is the Simulink model derived? In your report, you need to present (i) all the relevant equation derivations based on which you get the Simulink model, and (ii) the complete Simulink model.
(6 Marks)
c) For each block you have in your Simulink model, give details on all the settings, and the reasons for your settings.
(6 Marks)
(2) (5 Marks) Use the Simulink model to plot the response for the following input, f(t). The initial conditions are zero, i.e. 𝑥1 = 0, 𝑥̇1 = 0; 𝑥2 = 0, 𝑥̇2 = 0.
𝑡 0≤𝑡≤1 𝑓(𝑡) = {2 − 𝑡 1 < 𝑡 < 2
0 𝑡≥2
(II-A.2)
Hints:
Figure II-A.1. Two-mass suspension model.
a) The input to your model, can be implemented by using “Signal Builder” block.
b) To complete this question, you can also read and understand contents in my lecture
notes: Matlab Tutorial 3-2018-19, specially from pp. 21.
12 / 16
Question II-B, The motion of a single-mass suspension model of the vehicle suspension system.
(20 Marks)
Consider the single-mass suspension model shown in Figure II-B.1, which is typical representation of the vehicle suspension system. The spring and damper forces 𝑓 and
𝑓 have the nonlinear models. When the vehicle strikes a bump, the road surface profile 𝑑
can be represented by the trapezoidal function y(t) shown in Figure II-B.2. Now, the system model (derived from Newton’s law) becomes:
m𝑥̈ = 𝑓 (𝑦 − 𝑥) + 𝑓 (𝑦̇ − 𝑥̇) 𝑠𝑑
(II-B.1) where m = 400 kg, 𝑓 (𝑦 − 𝑥) is the is the nonlinear spring function shown in Figure
𝑠
II-B.3, and 𝑓 (𝑦̇ − 𝑥̇ ) is the nonlinear damper function given by: 𝑑
where 𝑣=𝑦̇−𝑥̇.
𝑓 (𝑣) = {−500|𝑣|1.2
𝑣 ≤ 0 𝑣 > 0
𝑑
50|𝑣|1.2
(1) (15 Marks) Develop a Simulink model of this system. In you report, you need to state clearly:
a) Answer a question: whether to use a state-variable representation, or a transfer- function representation, or nonlinear treatment, to create the corresponding model in Simulink, and why?
(3 Marks)
b) How is the Simulink model derived? In your report, you need to present (i) all the relevant equation derivations based on which you get the Simulink model, and (ii) the complete Simulink model.
(6 Marks)
c) For each block you have in your Simulink model, give details on all the settings, and the reasons for your settings.
(6 Marks)
(2) (5 Marks) Use the Simulink model to plot the response of this system, i.e. x and
(y-x) versus t. The initial conditions are zero, i.e. 𝑥 = 0 and 𝑥̇ = 0.
a) Plot the responses in Simulink using Scope. In your report, present the screen
shots.
(2 Marks)
b) Export the data into MATLAB workspace which enable you to use “plot” command to plot the responses. In your report, state (i) how do you export Simulink results to workspace, and, (ii) present the plots you have done using “plot” command, and show clearly the axis labels (i.e. for x and y axis, what quantities do you plot and units) and legends if applicable.
(3 Marks)
13 / 16
𝑠
(II-B.2)
Hints:
a) The road surface profile, y(t) as shown in Figure II-B.2, which is the input to your model, can be implemented by using “Signal Builder” block.
b) The spring function, as shown in Figure II-B.3, can be implemented by using “Look-Up Table” block.
c) The damper function, given by Equation (II-B.2), can be implemented by using “MATLAB function” block.
d) To complete this question, you can also read and understand contents in:
My lecture notes: Matlab Tutorial 3-2018-19, specially from pp. 47.
Textbook: Palm William, Introduction to MATLAB for Engineers (3rd Edn.), McGraw Hill Professional, Section 10.9, Chapter 10.
Figure II-B.1. Single-mass model of a vehicle suspension.
Figure II-B.2. Road surface profile. 14 / 16
Figure II-B.3. Nonlinear spring function.
15 / 16
Contribution Form – MATLAB Group Assignment Group____:
Name and Student ID number
Overall contribution (%)
Problem analysis and equation derivations, etc. (%)
MATLAB programming (%)
Report writing (%)
Comments
A 12345678
37
40
60
10
For example: Lead in the completion of the assignment and contribute the most in MATLAB programming
B 12345678
30
20
15
55
For example: Work closely with other group members, contributes particularly in Component 3.
C 12345678
28
35
20
30
For example: Work closely with other group members, contributes almost evenly in Component 1 and 3.
D 12345678
5
5
5
5
For example: He/she contributes very little.
Total
100%
100%
100%
100%
⚫ The parts highlighted in red are which you need to fill in.
⚫ The contribution form will be used to monitor the workload distribution within the group, and will also provide a basis for differentiate your marks (if there are substantial workload differences between team members).
16 / 16