EE4323 Industrial Control Systems
Homework Assignment 2
Nonlinear DC motor
Paulo R. Loma Marconi
June 13, 2017
1 Introduction
The objective is to simulate a nonlinear electro-mechanical system with thermal model and static
Coulomb friction. We use three ODE solvers, the embedded Matlab solver ode45, and two external
solvers, the 4th and 5th order Runge-Kutta algorithm ode45m, and the basic Euler algorithm eufix1.
2 Nonlinear model
The dynamics of the DC motor has two nonlinear parameters is,
R
A
i
A
+ L
A
˙
i
A
+ αω
1
= e
i
(t) (1)
J
1
˙ω
1
+ B
1
ω
1
− r
1
f
c
= αi
A
(2)
J
2
˙ω
2
+ B
2
ω
2
+ B
2C
sign(ω
2
) + r
2
f
c
= −τ
L
(3)
C
T M
˙
θ
M
+
(θ
M
− θ
A
)
R
T M
= i
2
A
R
A
(4)
where R
A
is the stationary resistance, L
A
is the stationary inductance, i
A
is the input stationary
current, α is the internal parameters, ω is the angular speed, e
i
(t) is the applied armature voltage,
B is the rotational viscous-damping coefficient, J is the moment of inertia, f
c
is the contact force
between two gears, r is the gear radius, and B
2C
is the static Coulomb friction. The thermal model
is similar to an electrical capacitor-resistor model with thermal capacity C
T M
, R
T M
is the resistive
losses to ambient temperature, θ
M
is the motor temperature, and θ
A
is the ambient temperature.
Now let us define the sate-vector differential equations: state vector x = [i
A
ω
2
θ
M
]
T
, and input
vector u = [e
i
τ
L
θ
A
]
T
.
For ω
1
= Nω
2
and N =
r
2
r
1
, eliminating f
c
we have
˙
i
A
= −
R
A
L
A
i
A
−
Nα
L
A
ω
2
+
1
L
A
e
i
(5)
˙ω
2
=
Nα
J
eq
i
A
−
B
eq
J
eq
ω
2
−
B
2C
J
eq
sign(ω
2
) −
1
J
e
q
τ
L
(6)
˙
θ
M
=
R
A
C
T M
i
2
A
−
1
C
T M
R
T M
θ
M
+
1
C
T M
R
T M
θ
A
(7)
where J
eq
= J
2
+ N
2
J
1
and B
eq
= B
2
+ N
2
B
1
.
1
For simulation purpose only we can simplify as
˙
i
A
= −a i
A
− b ω
2
+
1
L
A
e
i
(8)
˙ω
2
= c i
A
− d ω
2
− e sign(ω
2
) −
1
J
e
q
τ
L
(9)
˙
θ
M
= f i
2
A
− g θ
M
+ g θ
A
(10)
where a =
R
A
L
A
, b =
Nα
L
A
, c =
Nα
J
eq
, d =
B
eq
J
eq
, e =
B
2C
J
eq
, f =
R
A
C
T M
, and g =
1
C
T M
R
T M
.
3 Matlab Scripts
3.1 ODE solver ode45m
Listing 1 : ode45m
1 fu n c t i o n [ tout , yout ] = o d e 45m ( ypfun , t0 , tfinal , y0 , tol , trace )
2 % ODE45 Solve d ifferent i a l equation s , h i gher order meth o d .
3 % O DE45 in t e g rates a system of ordinary d i f ferenti a l equations using
4 % 4 th and 5 th order Runge - K utta formul a s .
5 % [T , Y ] = ODE45 (’ yprime ’ , T0 , Tfinal , Y0 ) integr a t e s the sy s t em of
6 % ordinary d i f ferenti a l equations de s cr ib ed by the M - file YPRIM E .M ,
7 % over the interval T0 to Tfinal , with initi a l conditio n s Y0 .
8 % [T , Y ] = ODE4 5 (F , T0 , Tfinal , Y0 , TOL , 1) uses tolerance TOL
9 % and d i s p l ay s st a t us while the integ r a tion proce e d s .
10 %
11 % I NPUT :
12 % F - S t r ing c ontaining name of user - supplied p r o b l e m descript i o n .
13 % Call : yprime = fun (t ,y ) w here F = ’fun ’.
14 % t - Time ( scalar ).
15 % y - So l u ti o n column - v ector .
16 % yprime - Returned d er ivative column - vector ; yprime (i ) = dy (i) / dt .
17 % t0 - In i t i a l v alue of t.
18 % tfinal - Final v alue of t.
19 % y0 - In i t i a l v alue column - vector .
20 % tol - The d e s ired ac c u r a c y . ( Default : tol = 1. e -6) .
21 % t race - If nonzero , each step is printed . ( Default : trace = 0) .
22 %
23 % OU T PUT :
24 % T - Re t u rn e d integr a t ion time points ( column - vector ) .
25 % Y - Re t u rn e d solut ion , one solution column - vector per tout - value .
26 %
27 % The r esult can be displayed by : plot ( tout , yout ) .
28 %
29 % See also ODE23 , ODEDEMO .
30
31 % C.B . Moler , 3 -25 -87 , 8 -26 -91 , 9 -08 -92.
32 % Copyright ( c) 1984 -94 by The MathWork s , Inc .
33
34 % The F e h l b e rg coeff i c ients :
35 alpha = [1/4 3/8 12/13 1 1/2] ’;
36 beta = [ [ 1 0 0 0 0 0]/4
37 [ 3 9 0 0 0 0]/32
38 [ 1932 -7200 7296 0 0 0]/2197
39 [ 8341 -32832 29440 -845 0 0]/4104
40 [ -6080 41040 -28352 9295 -5643 0 ] / 2 0 5 2 0 ] ’;
41 gamma = [ [902880 0 3953664 3 8 5 5 7 3 5 -1371249 27702 0 ] /7618 0 5 0
42 [ -2090 0 22528 21970 -15048 -27360]/752400 ] ’;
43 pow = 1/5;
44 if narg i n < 5, tol = 1. e -6; end
45 if narg i n < 6, trace = 0; end
46
47 % Init i a lizati o n
48 hmax = ( tfina l - t0 ) /16 ;
2
49 h = hmax /8;
50 t = t0 ;
51 y = y0 (:) ;
52 f = ze ros ( l e n gth ( y) ,6) ;
53 chunk = 128;
54 tout = z eros ( chunk ,1) ;
55 yout = z eros ( chunk , length (y ));
56 k = 1;
57 tout ( k) = t;
58 yout ( k ,:) = y . ’;
59
60 if tr ace
61 clc , t , h , y
62 end
63
64 % The main loop
65
66 while (t < t f inal ) & (t + h > t )
67 if t + h > tfinal , h = tfin a l - t; end
68
69 % Compu t e the s l o pes
70 temp = fe val ( ypfun , t , y );
71 f (: ,1) = temp (:) ;
72 for j = 1:5
73 temp = feval ( ypfun , t + al pha ( j)*h , y + h* f * beta (: , j) ) ;
74 f (: , j +1) = temp (:) ;
75 end
76
77 % Estimate the error and the ac c e p t able error
78 delta = norm ( h*f * ga mma (: ,2) ,’ inf ’);
79 tau = tol * max ( n orm (y , ’ inf ’ ) ,1.0) ;
80
81 % Up d a te the s o l u t i o n o nly if the error is a c c e pt able
82 if delt a <= tau
83 t = t + h;
84 y = y + h* f * gam ma (: ,1) ;
85 k = k +1;
86 if k > le n g th ( tout )
87 tout = [ tout ; zero s ( chunk ,1) ];
88 yout = [ yout ; zero s ( chunk , length (y )) ];
89 end
90 tout ( k ) = t ;
91 yout (k ,:) = y . ’;
92 end
93 if trac e
94 home , t , h , y
95 end
96
97 % Up d a te the step size
98 if delt a ~= 0.0
99 h = min ( hmax , 0.8* h *( tau / delta )^ pow );
100 end
101 end
102
103 if (t < t f inal )
104 disp ( ’ S i ngularit y li kely . ’)
105 t
106 end
107
108 tout = tout (1: k) ;
109 yout = yout (1: k ,:) ;
3.2 ODE solver eufix1
Listing 2 : eufix1
3
1 f u n c t i o n [ tout , xout ] = e u f ix1 ( dxfun , tspan , x0 , stp , tra ce )
2 % EUFIX1 Sol ve o r d i na r y state - vector di f f e rential e quations , low orde r m ethod .
3 % EU F IX1 i ntegrates a set of ODEs xdot = f (x , t ) usin g the most
4 % elementa r y Euler algorithm , witho u t step - size control .
5 %
6 % CALL :
7 % [t , x] = eufix1 ( ’ dxfun ’, tspan , x0 , stp , trace )
8 %
9 % I NPUT :
10 % d xfun - String con t a i ning name of user - supplied p r o b l e m descrip t i on .
11 % Call : xdot = m odel (t ,x ) c oded in fname . m => dxfu n = ’ fname ’.
12 % t - T ime ( scalar ) .
13 % x - Solution column - vector at time t .
14 % xdot - Returned d e r ivative column - vector ; xdot = dx / dt .
15 % t span - Range of t for the desired solution ; tspan = [ t0 tf ].
16 % tf - F inal valu e of t.
17 % x0 - In i t i a l v alue column - vector .
18 % stp - The specified i ntegrati o n step ( default : stp = 1. e -2) .
19 % t race - If nonzero , each step is printed ( default : trace = 0) .
20 %
21 % OU T PUT :
22 % t - Re t u rn e d integr a t ion time points ( row - vecto r ).
23 % x - Re t u rn e d solut ion , one column - vector per tout - value .
24 %
25 % Disp l a y re s ult by : plot (t , x ) or plot (t , x (: ,2) ) or pl ot (t , x (: ,2) , x (: ,5) ) .
26
27 % Init i a lizati o n
28 if narg i n < 4, stp = 1. e -2; disp (’H = 0.02 by default ’ ) ; end
29 if narg i n < 5, trace = 0; end %% d isable trace if not r e q u e st ed
30 t0 = tsp an (1) ; tf = ts pan (2) ;
31 if tf < t0 , error ( ’tf < t0 ! ’); return ; end % % chec k for glaring error
32 t = t0 ;
33 h = stp ;
34 x = x0 (:) ;
35 k = 1;
36 tout ( k) = t; % initializ e ou tput arrays
37 xout ( k ,:) = x . ’;
38 if tr ace
39 clc , t , h , x
40 end
41
42 % The main loop
43
44 while (t < tf )
45 if t + h > tf , h = tf - t; end
46 % Compu t e the derivati v e
47 dx = feval (dxfun , t , x ) ; dx = dx (:) ;
48 % Up d a te the s o l u t i o n ( wi th no c heck on e rror )
49 t = t + h ;
50 x = x + h * dx ;
51 k = k +1;
52 tout (k ) = t;
53 xout (k ,:) = x . ’;
54 if trac e
55 home , t , h , x , dx
56 end
57 end
58 if (t < tf ) % if true , somethin g bad h a p p en e d !
59 disp ( ’ S i ngularit y or mod e l i n g error likely . ’)
60 t
61 end
62 % ... here is the outpu t ( tout in row vector form )
63 tout = tout (1: k) ;
64 xout = xout (1: k ,:) ;
4
3.3 Nonlinear model
In line 37 and 38, the input e
i
can be changed from constant input to sinusoidal input.
Listing 3 : Nonlinear model
1 f u n c t i o n xdot = ass t 0 2 _2017 (t ,x )
2 global E_0 Tau_L0 T_Amb B_2C
3
4 % m otor parameters , Nachtiga l , Table 16.5 p . 663
5
6 J_1 = 0.0035; % in * oz *s ^2/ rad
7 B_1 = 0.064; % in * oz *s / rad
8
9 % electric a l / mechanica l relations
10 K_E = 0.1785; % back emf co efficient , e_m = K_E * omega_m ( K_E = alpha * omega )
11 K_T = 141.6* K_E ; % torqu e co e f f i c ., in English units K_T is not = K_E ! ( K_T = a lpha *
iA )
12 R_A = 8.4; % Ohms
13 L_A = 0.0084; % H
14
15 % gear - train and load para m e t ers
16 J_2 = 0.035; % in * oz *s ^2/ rad % 10 x motor J
17 B_2 = 2.64; % in * oz *s / rad ( viscous )
18 N = 8; % mot or / l oad gear rati o ; omega_1 = N omega_2
19
20 % Ther m a l model pa r a m et ers
21 R_TM = 2.2; % K elvin / Watt
22 C_TM = 9/ R_TM ; % Watt - sec / Kelvi n (-> 9 sec time constant - fast !)
23
24 Jeq = J_2 + N *2* J_1 ;
25 Beq = B_2 + N ^2* B_1 ;
26 a = R_A / L_A ;
27 b = K_E *N / L_A ;
28 c = N* K_T / Jeq ;
29 d = Beq / Jeq ;
30 e = B_2C / Jeq ;
31 f = R_A / C_TM ;
32 g = 1/( C_TM * R_TM );
33
34 if t < 0.05
35 e_i = 0;
36 else
37 % e_i = E_0 ;
38 e_i = E_0 * sin (5 *(2* pi ) *( t - 0.05) ) ;
39 end
40 if t < 0.2
41 Tau_L = 0;
42 else
43 Tau_L = Tau _ L0 ;
44 end
45
46 xdot (1) = - a *x (1) - b *x (2) + e_i / L_A ;
47 xdot (2) = c * x (1) -d * x (2) -e * sign (x (2) ) - Tau _L / Jeq ;
48 xdot (3) = f * x (1) ^2 - g*x (3) + g* T_Amb ;
49 xdot = xdot (:) ; % force column vector
3.4 Main
Change the input values in line 5-8, the input_type E
0
to constant or sinusoidal, and the step size in
line 9.
Note: This script is an example only. In the Simulation Results section we will analyze different
scenarios.
Listing 4 : Main
5
1 clear variables ; close all ; clc ;
2 global E _0 Tau _L 0 T_ Amb B _2 C;
3
4 E _0 = 120; % [V] 120
5 Tau _ L 0 = 80; % [ N .m ] 80
6 T _ Amb = 18; % [ deg ] 18
7 B _2 C = 300; % [ N ] 80/300
8
9 t 0 = 0; t f i nal = 0.3; step = 1e -4;
10 x 0 = [0; 0; 0]; % i nitial co nditions
11
12 input _ type = 0; % 0= constant , 1= si n u s o idal
13 % % ode 45 vs ode 45 m vs eu fix 1
14
15 timer = clo ck ;
16 [ t1 , x 1] = ode 45( ’ a sst 02_2017 ’ ,[ t0 , tfi n a l ], x 0) ;
17 % [t1 , x 1] = ode 45 m (’ asst 02_2017 ’ , t0 , tfinal , x 0, step ) ;
18 Tsim 1 = etime ( clock , tim er ) ; % integr a t i on ti me
19 Len 1 = le n g th ( t 1) ; % number of time - steps
20
21 timer = clo ck ;
22 [ t2 , x 2] = ode 45 m ( ’ as st 02_2017 ’ , t0 , tfinal , x 0, step ) ;
23 Tsim 2 = etime ( clock , tim er ) ; % integr a t i on ti me
24 Len 2 = le n g th ( t 2) ; % number of time - steps
25
26 timer = clo ck ;
27 [ t3 , x 3] = eufix 1( ’ asst 02_2017 ’ ,[ t 0 tfi n a l ], x0, step );
28 Tsim 3 = etime ( clock , tim er ) ; % integr a t i on ti me
29 Len 3 = le n g th ( t 3) ; % number of time - steps
30
31 % % Relative error
32
33 % relative error at max current : ode 45 vs eufix 1
34 max _ iA _ ode 45 = max ( x 1(: ,1) );
35 max _ iA _ eufix 1 = max ( x 3(: ,1) ) ;
36 max _ iA _ error = 100* abs ( ( max _ iA _ ode 45 - max _ iA _ eufix 1) / max _ iA _ ode 45 ) ;
37
38 % relative error at max angular velocity : ode 45 vs e ufix 1
39 max _ omega 2_ ode 45 = max (x 1(: ,2) ) ;
40 max _ omega 2_ e ufix 1 = max ( x 3(: ,2) );
41 max _ omega 2_ e rror = 100* abs ( ( max _ omega 2_ ode 45- max _ omeg a 2_ e ufix 1)/ max _ o mega 2_ ode 45 );
42
43 % % Plotting
44 if in put _ type == 0
45 %% Constant inpu t e _i = E0
46 figure ;
47 subplot (3 ,1 ,1) ;
48 plot (t 1,x 1(: ,1) ,t 2 ,x 2(: ,1) ,’- - ’ ,t3 , x 3(: ,1) ,’-. ’ ,’ LineWidth ’ ,1.5) ;
49 ti tle ([ ’ Non l i n e a r DC motor with thermal model ,
$
B _{2 C }=
$
’ , num 2 str ( B _2 C)] , ’
Int e rpreter ’, ’ Latex ’) ;
50 yla b e l (’
$
i_ A
$
[A ] ’, ’ I n terpreter ’,’ Latex ’) ;
51 leg e n d ([ ’ ode 45: ’, num 2 str ( Tsim 1) ,’ [ s ] ’] ,[ ’ ode 45 m : ’, num 2 str ( Tsim 2) ,’ [ s ] ’] ,[ ’
euf ix 1: ’, num 2 str ( Tsim 3) ,’ [s ] ’]) ;
52 grid on ;
53
54 subplot (3 ,1 ,2) ;
55 plot (t 1,x 1(: ,2) ,t 2 ,x 2(: ,2) ,’- - ’ ,t3 , x 3(: ,2) ,’:’ , ’ LineWidth ’ ,1.5) ;
56 yla b e l ( ’
$
\ omega _2
$
[ rad /s ] ’, ’ Interpr eter ’ , ’ Latex ’) ;
57 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
58 grid on ;
59
60 subplot (3 ,1 ,3) ;
61 plot (t 1,x 1(: ,3) ,t 2 ,x 2(: ,3) ,’- - ’ ,t3 , x 3(: ,3) ,’:’ , ’ LineWidth ’ ,1.5) ;
62 xla b e l (’ Time [s ] ’,’ I n terpreter ’, ’ Latex ’) ;
63 yla b e l ( ’
$
\ theta _ M
$
[ deg ]’ , ’ Interp r eter ’ , ’ Latex ’) ;
64 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
65 grid on ;
66
6
67 % print ( ’../ asst 02_2017/ E 0_ ode 45 - ode 45m - euf ix 1_1 e -3. png ’, ’- dpng ’ , ’-r 300 ’) ; % Save
as PNG with 300 DPI
68
69 figure ;
70 subplot (2 ,1 ,1) ;
71 plot (t 1,x 1(: ,1) ,t 2 ,x 2(: ,1) ,’- - ’ ,t3 , x 3(: ,1) ,’-. ’ ,’ LineWidth ’ ,1.5) ;
72 ti tle ([ ’ Non l i n e a r DC motor with thermal model ,
$
B _{2 C }=
$
’ , num 2 str ( B _2 C)] , ’
Int e rpreter ’, ’ Latex ’) ;
73 yla b e l (’
$
i_ A
$
[A ] ’, ’ I n terpreter ’,’ Latex ’) ;
74 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’) ;
75 axis ([0.05 0.0 7 - inf inf ]) ;
76 text (0.058 ,5.5 ,[ ’ Rela t i v e error at max
$
i _{ A }
$
= ’ , num 2 str ( max _ iA _ error ) ,’
$
\ %
$
’] , ’ Interprete r ’ ,’ Latex ’) ;
77 grid on ;
78
79 subplot (2 ,1 ,2) ;
80 plot (t 1,x 1(: ,2) ,t 2 ,x 2(: ,2) ,’- - ’ ,t3 , x 3(: ,2) ,’:’ , ’ LineWidth ’ ,1.5) ;
81 yla b e l ( ’
$
\ omega _2
$
[ rad /s ] ’, ’ Interpr eter ’ , ’ Latex ’) ;
82 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
83 axis ([0.05 0.0 7 - inf inf ]) ;
84 text (0.058 ,40 ,[ ’ R e l a t i v e error at max
$
\ omega _{2}
$
= ’ , num 2 str ( max _ o mega 2_ error )
,’
$
\ %
$
’] , ’ Interprete r ’ ,’ Latex ’) ;
85 grid on ;
86
87 % print ( ’../ asst 02_2017/ E 0_ ode 45 - ode 45m - euf ix 1_1 e -3_ zoom . png ’, ’ - dpng ’,’ -r 300 ’) ; %
Save as PNG with 300 DPI
88
89 figure ;
90 plot (t 1,x 1(: ,3) ,t 2 ,x 2(: ,3) ,’- - ’ ,t3 , x 3(: ,3) ,’:’ , ’ LineWidth ’ ,1.5) ;
91 ti tle ( ’ Mo tor t e mperature
$
\ theta _ M
$
over
$
80[ s ]
$
’ , ’ Interpret e r ’ ,’ Latex ’) ;
92 xla b e l (’ Time [s ] ’,’ I n terpreter ’, ’ Latex ’) ;
93 yla b e l ( ’
$
\ theta _ M
$
[ deg ]’ , ’ Interp r eter ’ , ’ Latex ’) ;
94 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
95 grid on ;
96
97 % print ( ’../ asst 02_2017/ thetaM _ode 45- ode 45 m - eu fix 1_1e -3. png ’, ’- dpng ’, ’ -r 300 ’) ; %
Save as PNG with 300 DPI
98
99 elseif input _type == 1
100 %% Sinusoida l input e_ i
101
102 figure ;
103 subplot (3 ,1 ,1) ;
104 plot (t 1,x 1(: ,1) ,t 2 ,x 2(: ,1) ,’- - ’ ,t3 , x 3(: ,1) ,’-. ’ ,’ LineWidth ’ ,1.5) ;
105 ti tle ([ ’ Non l i n e ar DC motor with thermal model ,
$
B _{2 C }=
$
’ , num 2 str ( B _2 C)] , ’
Int e rpreter ’, ’ Latex ’) ;
106 yla b e l (’
$
i_ A
$
[A ] ’, ’ I n terpreter ’,’ Latex ’) ;
107 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
108 grid on ;
109
110 subplot (3 ,1 ,2) ;
111 plot (t 1,x 1(: ,2) ,t 2 ,x 2(: ,2) ,’- - ’ ,t3 , x 3(: ,2) ,’:’ , ’ LineWidth ’ ,1.5) ;
112 yla b e l ( ’
$
\ omega _2
$
[ rad /s ] ’, ’ Interpr eter ’ , ’ Latex ’) ;
113 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
114 grid on ;
115
116 subplot (3 ,1 ,3) ;
117 plot (t 1,x 1(: ,3) ,t 2 ,x 2(: ,3) ,’- - ’ ,t3 , x 3(: ,3) ,’:’ , ’ LineWidth ’ ,1.5) ;
118 xla b e l (’ Time [s ] ’,’ I n terpreter ’, ’ Latex ’) ;
119 yla b e l ( ’
$
\ theta _ M
$
[ deg ]’ , ’ Interp r eter ’ , ’ Latex ’) ;
120 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
121 grid on ;
122
123 % print ( ’../ asst 02_2017/ sinE 0_ ode 45 - ode 45m - euf ix 1_1 e -4. png ’, ’ - dpng ’, ’-r 300 ’) ; %
Save as PNG with 300 DPI
124
125 figure ;
126 subplot (3 ,1 ,1) ;
7
127 plot (t 1,x 1(: ,2) ,t 2 ,x 2(: ,2) ,’- - ’ ,t3 , x 3(: ,2) ,’-. ’ ,’ LineWidth ’ ,1.5) ;
128 ti tle ([ ’ Stiction behaviour on
$
\ omega _2
$
,
$
B _{2 C }=
$
’ , num 2 str ( B _2 C)] , ’
Int e rpreter ’, ’ Latex ’) ;
129 yla b e l ( ’
$
\ omega _2
$
[ rad /s ] ’, ’ Interpr eter ’ , ’ Latex ’) ;
130 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
131 axis ([ 0 . 1 48 0.157 -0.6 0.4] ) ;
132 grid on ;
133
134 subplot (3 ,1 ,2) ;
135 plot (t 1,x 1(: ,2) ,t 2 ,x 2(: ,2) ,’- - ’ ,t3 , x 3(: ,2) ,’:’ , ’ LineWidth ’ ,1.5) ;
136 yla b e l ( ’
$
\ omega _2
$
[ rad /s ] ’, ’ Interpr eter ’ , ’ Latex ’) ;
137 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
138 axis ([ 0 . 1 48 0.157 -0.015 0.010]) ;
139 grid on ;
140
141 subplot (3 ,1 ,3) ;
142 plot (t 1,x 1(: ,2) ,t 2 ,x 2(: ,2) ,’- - ’ ,t3 , x 3(: ,2) ,’-. ’ ,’ LineWidth ’ ,1.5) ;
143 xla b e l (’ Time [s ] ’,’ I n terpreter ’, ’ Latex ’) ;
144 yla b e l ( ’
$
\ omega _2
$
[ rad /s ] ’, ’ Interpr eter ’ , ’ Latex ’) ;
145 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
146 axis ([ 0 . 1 48 0.157 -11 e -5 5e -5]) ;
147 grid on ;
148
149 % print ( ’../ asst 02_2017/ sinE 0_ ode 45 - ode 45m - euf ix 1_1 e -4_ zoom . png ’, ’- dpng ’ , ’- r
300 ’) ; % Save as PNG with 300 DPI
150
151 figure ;
152 plot (t 1,x 1(: ,3) ,t 2 ,x 2(: ,3) ,’- - ’ ,t3 , x 3(: ,3) ,’:’ , ’ LineWidth ’ ,1.5) ;
153 ti tle ( ’ Mo tor t e mperature
$
\ theta _ M
$
over
$
80[ s ]
$
’ , ’ Interpret e r ’ ,’ Latex ’) ;
154 xla b e l (’ Time [s ] ’,’ I n terpreter ’, ’ Latex ’) ;
155 yla b e l ( ’
$
\ theta _ M
$
[ deg ]’ , ’ Interp r eter ’ , ’ Latex ’) ;
156 leg e n d (’ ode 45’ , ’ ode 45m ’,’ e ufix 1 ’ ,’ Location ’ ,’ sout heast ’) ;
157 grid on ;
158
159 % print ( ’../ asst 02_2017/ sinE 0_ thetaM _ ode 45 - ode 45m - eufix 1_1 e -4. png ’,’ - dpng ’ ,’- r
300 ’) ; % Save as PNG with 300 DPI
160 end
4 Simulation scenarios
Two types of scenarios are simulated, Table 1. The first one is submitted to a constant input, low
stiction, and two step sizes. The second scenario is more interesting because we study the behaviour
due to a sinusoidal input which emulates the reversing mode of the motor at 5[Hz] with higher stiction.
Both scenarios have load torque at t = 0.2[s].
Scenario 1 Scenario 2
ode45m step size 1 × 10
−3
1 × 10
−4
1 × 10
−4
eufix1 step size 1 × 10
−3
1 × 10
−4
1 × 10
−4
ode45 step size auto auto
e
i
E
0
E
0
sin[5(2π)(t − 0.05)]
E
0
120 [V ] 120 [V ]
τ
L
80 [Nm] at t = 0.2[s] 80 [N m] at t = 0.2[s]
θ
A
18 [
°
C] 18 [
°
C]
B
2C
80 [N] 300 [N]
Table 1: Scenario 1 and 2.
8
5 Simulation Results
Scenario 1
The result in Fig. 1 shows the output states due to constant input E
0
= 120, and B
2C
= 80. The
current overshoot at 0.05[s] is due to the inertia that the motor has to overcome. After the inertia is
broken, the current i
A
drops down to a constant value. The load torque τ
L
is applied at 0.2[s] which
produces the increment in the current and the decrement in the angular velocity. Also, eufix1 solves
the system with noticeable error, this result is analyzed later.
Figure 1: Scenario 1: step size 1 × 10
−3
The time simulation of each solver indicates that ode45 is 10 times slower than ode45m and 30
times slower than eufix1.
ode45 ode45m eufix1
simulation time [s] 0.267 0.025 0.009
number of time steps 71769 15133 80001
Table 2: Scenario 1: simulation time and number of steps.
The temperature of the motor θ
M
increases linearly and reaches steady-state at 45[s] approximately,
which indicates that the motor won’t reach unsafe temperatures.
9
Figure 2: Scenario 1: θ
M
over 80[s]
Although, eufix1 is the fastest solver, with step size of 1 × 10
−3
, eufix1 outputs the worst
performance. The result can be improved if the steps size is decreased to 1 × 10
−4
. Fig. 3 and Fig. 4
show the relative error at max current and max angular velocity between ode45 and eufix1.
Figure 3: Scenario 1: step size 1 × 10
−3
10
Figure 4: Scenario 1: step size 1 × 10
−4
Scenario 2
In this scenario we submitted the DC motor to high stiction B
2C
= 300 and sinsusoidal input at 5[Hz]
which simulates the reversing mode. Table 3 shows the output for the three solvers showing that ode45
is the slowest by far.
ode45 ode45m eufix1
simulation time [s] 517.798 3.463 0.093
number of time steps 13559913 9647 3002
Table 3: Scenario 2: simulation time and number of steps.
Fig. 5 shows the simulation output. The relevant result is the behavior of the system around the
(nonlinear) stiction. ω
2
sticks at 0.15[s] and 0.25[s] due to B
2C
.
11
Figure 5: Scenario 2: reversing mode.
Even though the solvers were able to solve the dynamics with stiction, ode45 took too much time
to overcome this nonlinearity. Fig. 6 shows the stiction with three different zoom levels for each solver.
The fastest but with more integration step error is eufix1. ode45m has less error ±0.01, and finally
ode45 solves with the minimum error, around ±10 × 10
−5
. In conclusion, ode45m is the best choice
against the rest because it can obtain the solution with low error and with decent speed.
12
Figure 6: Scenario 2: stiction behavior at 0.148 ≤ t ≤ 0.157.
The motor temperature θ
M
in reversing mode reaches the steady-state at 45[s] approximately which
indicates the motor won’t suffer overheat due to stiction.
Figure 7: Scenario 2: motor temperature.
13
