Newton’s Cradle Example

Consider a simplified model of the Newton’s cradle consisting of a pair of pendulums with mass m_L, m_R as shown in Figure 1:

figure1

The Hybrid Model for the system is given by

f(x) =\left[  \begin{array}{c}  x_2\\  -\frac{\gamma}{l} \cdot \sin(x_1)-b_L \cdot x_2\\  x_4\\  -\frac{\gamma}{l} \cdot \sin(x_3)-b_R \cdot x_4\\  \end{array}  \right]

where \gamma is the gravity constant, b_L= \frac{\widetilde{b_L}}{m_L l^2}, b_R= \frac{\widetilde{b_R}}{m_R l^2}; m_L l^2=m_R l^2 are the Pendulum’s moment of inertia and \widetilde{b_L}, \widetilde{b_R} are the coefficients of viscous friction torque for each pendulum.

C:=\{x: x_1<x_3\} \cup \{x: x_1=x_3,x_2<x_4\}

g(x) =\left[  \begin{array}{c}  x_1\\  (\lambda -(1-\lambda)e) x_2+( 1- \lambda) (1+e) x_4 \\  x_3\\  \lambda (1+e) x_2 + (1-\lambda-\lambda e) x_4 \\  \end{array}  \right]

where \lambda=\frac{m_L}{m_L+m_R} is the mass ratio and e \in (0,1) is the restitution coefficient.

D:=\{x: x_1=x_3, x_2>x_4\}

Source code for Simulation: script to run the simulation (run.m)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab M-file       Project: HyEQ Toolbox  @ Hybrid Dynamics and Control
% Lab, http://www.u.arizona.edu/~sricardo/index.php?n=Main.Software
%
% Filename: run
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function run
%constants
b_Ltilde = 0.25;
b_Mtilde = 0.25;
global gam l b_L b_R m_L m_R e;

gam = 9.81;
l =  1;
m_L = 1;
m_R = 2;
b_L = b_Ltilde/m_L*l^2;
b_R = b_Ltilde/m_R*l^2;
e = 0.75;

% initial conditions
%x = [x1; x2; x3; x4]
% x0 = [-pi/4; 0; pi/4; 0]; % case1
% x0 = [0;0; pi/4; 0]; % case2
x0 = [pi/4;0; pi/2; 0]; % case3

% simulation horizon
TSPAN=[0 10];
JSPAN = [0 30];

% rule for jumps
% rule = 1 -> priority for jumps
% rule = 2 -> priority for flows
rule = 1;

options = odeset('RelTol',1e-4,'MaxStep',1e-3);

% simulate
[t j x] = HyEQsolver( @f,@g,@C,@D,...
    x0,TSPAN,JSPAN,rule,options);

% state
x1 = x(:,1);
x2 = x(:,2);
x3 = x(:,3);
x4 = x(:,4);

% plot solution
figure(1)
clf
% hold on
subplot(2,1,1),plotflows(t,j,x1)
grid on
ylabel('x1','rot',0)
xlabel('t')
axis([0 10 -1.75 1.75])
set(gca,'ytick',-pi/2:pi/4:pi/2)
set(gca,'YTickLabel',{'-pi2','rem','0','rem','pi2'})

subplot(2,1,2),plotflows(t,j,x2)
grid on
ylabel('x2','rot',0)
xlabel('t')
axis([0 10 -3.8 3.8])
axis([0 10 -4.5 4.5]) %case 1 and 2 and 3
set(gca,'ytick',-pi:pi/2:pi)
set(gca,'YTickLabel',{'-pi','rem','0','rem','pi'})

figure(2)
clf
% hold on
subplot(2,1,1),plotflows(t,j,x3)
grid on
ylabel('x3','rot',0)
xlabel('t')
% axis([0 10 -1 1])
axis([0 10 -1.75 1.75])
set(gca,'ytick',-pi/2:pi/4:pi/2)
set(gca,'YTickLabel',{'-pi2','rem','0','rem','pi2'})

subplot(2,1,2),plotflows(t,j,x4)
grid on
ylabel('x4','rot',0)
xlabel('t')
% axis([0 10 -2.5 2.5])
axis([0 10 -4.5 4.5]) %case 1 and 2 and 3
set(gca,'ytick',-pi:pi/2:pi)
set(gca,'YTickLabel',{'-pi','rem','0','rem','pi'})
end

% % plot hybrid arc
% plotHybridArc(t,j,x)
% xlabel('j')
% ylabel('t')
% zlabel('x')

Flow map (f.m)

function xdot = f(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab M-file
% Project: AME549_HW1_1.3
%
% Name: f.m
%
% Description: Flow map
%
% Version: 1.0
% Required files: -
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global gam l b_L b_R

% state
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);

% flow map
x1dot = x2;
x2dot = -gam/l*sin(x1) - b_L*x2;
x3dot = x4;
x4dot = -gam/l*sin(x3) - b_R*x4;

xdot = [x1dot; x2dot; x3dot; x4dot];
end

Flow set (C.m)

function [value discrete] = C(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab M-file       Project: HyEQ Toolbox  @ Hybrid Dynamics and Control
% Lab, http://www.u.arizona.edu/~sricardo/index.php?n=Main.Software
%
% Filename: C.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Description: Flow set
% Return C if outside of D, and 1 if inside C
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% state
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);

if (x1 >= x3) && (x2 < x4)
    value = 1;
elseif x1 < x3
    value = 1;
else
    value = 0;
end
end

Jump map (g.m)

function xplus = g(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab M-file       Project: HyEQ Toolbox  @ Hybrid Dynamics and Control
% Lab, http://www.u.arizona.edu/~sricardo/index.php?n=Main.Software
%
% Filename: g.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Project: Simulation of a hybrid system (bouncing ball)
% Description: Jump map
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global m_L m_R e

% state
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);

% constant
lam = m_L/(m_L + m_R);

% jump map
x1plus = x1;
x2plus = (lam - (1-lam)*e)*x2 + (1-lam)*(1+e)*x4;
x3plus = x3;
x4plus = lam*(1+e)*x2 + (1-lam - lam*e)*x4;

xplus = [x1plus; x2plus; x3plus; x4plus];

end

Jump set (D.m)

function xplus = g(x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab M-file       Project: HyEQ Toolbox  @ Hybrid Dynamics and Control
% Lab, http://www.u.arizona.edu/~sricardo/index.php?n=Main.Software
%
% Filename: g.m
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Project: Simulation of a hybrid system (bouncing ball)
% Description: Jump map
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global m_L m_R e

% state
x1 = x(1);
x2 = x(2);
x3 = x(3);
x4 = x(4);

% constant
lam = m_L/(m_L + m_R);

% jump map
x1plus = x1;
x2plus = (lam - (1-lam)*e)*x2 + (1-lam)*(1+e)*x4;
x3plus = x3;
x4plus = lam*(1+e)*x2 + (1-lam - lam*e)*x4;

xplus = [x1plus; x2plus; x3plus; x4plus];

end

Simulation Results are shown for \widetilde{b_L} =\widetilde{b_R}= 0.25, \gamma : 9.81 [m/s^2], l=1 [m], m_L = 1; m_R = 2; e = 0.75 and initial conditions: x_0=\begin{bmatrix}          \pi/4 & 0 & \pi/2 & 0    \\[0.3em]            \end{bmatrix}^T

sim1 sim2

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s