# How can I modeling a simple haptic interface system to obtain results similar to the paper “Time-Domain Passivity Control of Haptic Interfaces”?

I am trying to simulate the experiments in the paper "Time-Domain Passivity Control of Haptic Interfaces", Hanaford and Ryu, 2002, IEEE transactions on robotics and automation vol .18, No1 about the simulation 2, in Fig.8, and its results are fig.9 and fig.10.

I don't know an exact model of fig.9 which I can draw again in Matlab/simulink to run. I'm having trouble visualizing the mathematical model explaining the input and output of the system so that I can put it in a MATLAB code and perform the simulation. I am using constants like the image.

How can I do to get similar results? When I try to emulate something I get inconsistent results and I have trouble putting the mathematical model in lines of code.

I tried drawing spring and damper as in fig.6 of the same paper, but I don't know which is my input (force, velocity or position) and also the type of signal of input (sine, step,...) to get a similar result as the paper, or how I can insert the initial conditions on the code editor or on the Simulink.

I added 3 images, one is the parameters that I chose and another is the model which I drew.

I am really confused about initial condition of position.

In here, for the Virtual Environment: damping = -50 (I want to check the case without Passivity controller) k = 30000.

I am sending the simulation if it help. Thanks a lot :D

Matlab editor code: (sorry for the mess, I don't know how to put the modeling in the code) :/

%%%% Simulation Haptic Interface System %%%%

% - Impedance-based environment

clc, clear

%% Parameters

freq = 1000;           % Sampling rate simulation
delta_T = 1/freq;      % Sampling period simulation
t = 0:delta_T:10;      % Simulation time

%- Human Operator/ Disturbance (HO) -> mass + spring
M_HO = 0.1;            % Mass HO system (? kg)
k_HO = 10000;          % Spring HO constant  (? kN/m)

%- Haptic device (Interface) (HI)
M_HI = 0.1;            % Mass HI system (? kg)
k_HI = 150000;         % Spring HI constant  (? kN/m)
b_HI = 50;             % Damping HI constant (? Ns/m)

%- Delayed Virtual Environment (VE)
k_VE = 30000;          % Spring VE constant  (30 kN/m)
b_VE = - 50 ;          % Damping VE constant (-50 Ns/m)

freq_VE = 200/3;       % Sampling rate VE (66.67 Hz)
delta_Tve = 1/freq_VE; % VE system operates at a realativily slow samplig rate (15 ms)

%% Inicialization

%- Human Operator/ Disturbance (HO) -> mass + spring
x_HO = [0];            % Human position
dx_HO = [0];
ddx_HO = [0];
f_HO = [0];            % Force output (? N)

%- Haptic device (Interface) (HI)
x_HI = [0];            % Human disturbance input - position
dx_HI = [0];
ddx_HI = [0];
f_HI = [0];            % Force output (? N)

%- Series PC (alpha)
W = [0 0];             % Passivity Observer (We set up the PO to monitor only the
%                     virtual environment and the PC)
alpha = [0 0];         % Passivity Controller
f_PC = [0];

%- Delayed Virtual Environment (VE)
x_VE = [0];            % Haptic Interface disturbance input - position
dx_VE = [0];
ddx_VE = [0];
f_VE = [0];            % Force output (? N)

n = 2;

ve = [0];
fe = [0];
%% without PC (the system is highly unstable when driven to contact)

sum_ve = 0;
ve_inte = 0;

for t1 = 0.0000:delta_T:10

%- Human Operator
if t1 < 4
x_HO(:,end+1) = 0.140;                         % Human input (140 mm)
else
x_HO(:,end+1) = 0.090;                         % Human input (90 mm)
end

dx_HO(:,end+1) = (x_HO(n) - x_HO(n-1))/delta_T;
ddx_HO(:,end+1) = (dx_HO(n) - dx_HO(n-1))/delta_T;

f_HO(:,end+1) = M_HO*ddx_HO(n);

%- Haptic device
x_HI(:,end+1) = (f_HO(n) + k_HO*x_HO(n))/k_HO;

%     f_HI(:,end+1) = f_HO(n);                           % f_h

dx_HI(:,end+1) = (x_HI(n) - x_HI(n-1))/delta_T;    % v_h
ddx_HI(:,end+1) = (dx_HI(n) - dx_HI(n-1))/delta_T;

%- Virtual Environment

ve(:,end+1) = dx_HO(n) - dx_HI(n);

ve_inte = ve_inte + ve(n);

%     ve(n) = (x_VE(n) - x_VE(n-1))/delta_Tve;
x_VE(:,end+1) = ve_inte*delta_Tve;

dx_VE(:,end+1) = (x_VE(n) - x_VE(n-1))/delta_T;

fe(:,end+1) = k_VE*x_VE(n) + b_VE*ve(n);

n = n+1;
end

%  The signal flow should be as follows: haptic device
% --> Virtual Environment (f = ke + bv, where e = penetration depth inside
% the virtual environment/wall, v = haptic device velocity,
% k and b = stiffness and damping gains)
% --> PO/PC --> haptic device.

%% with added PC (the system achieves stable contact after about 3 bounces)

sum_ve = 0;
ve_inte = 0;

for t1 = 0.0000:delta_T:10

%- Human Operator
if t1 < 4
x_HO(:,end+1) = 0.140;                         % Human input (140 mm)
else
x_HO(:,end+1) = 0.090;                         % Human input (90 mm)
end

dx_HO(:,end+1) = (x_HO(n) - x_HO(n-1))/delta_T;
ddx_HO(:,end+1) = (dx_HO(n) - dx_HO(n-1))/delta_T;

f_HO(:,end+1) = M_HO*ddx_HO(n);

%- Haptic device

%- Virtual Environment

f_PC(:,end+1) = alpha(n)*ve(n);

f_VE(:,end+1) = fe(n) + alpha(n)*ve(n);

f_HI(:,end+1) = f_VE(n) + M_HI*ddx_HI(n) + b_HI*dx_HI(n) + k_HI*x_HI(n) ...
- b_HI*dx_VE(n) - k_HI*x_VE(n);

%- Series PC

W(:,end+1) = W(n-1) + (fe(n)*ve(n) + alpha(n-1)*(ve(n-1))^2);

sum_ve = sum_ve + ve(n)^2;

if W(n) < -b_HI*sum_ve
alpha(:,end+1) = (W(n) + b_HI*sum_ve)/(ve(n)^2);
else
alpha(:,end+1) = 0;
end

n = n+1;

end

%% Graphics


Best regards,

leonardo.fagundes@ufv.br

Edit1: simulation 2, in Fig.8, and its results are fig.9 and fig.10

edit2: In the first simulation, the response for simple virtual wall system when driven by a sinusoidal velocity profile generates energy when damping is negative (the system becomes active). When wall damping is still negative and passivity controller is operating, dissipation is constrained to be positive and the system is stable. I figured, in the second simulation, I should also use the same damping gain. The authors quote that in this simulation "The VE includes a spring constant of 30 kN/m and operates at a relatively slow sampling rate of 66.67 Hz (15 ms)"

Therefore, damping and spring constant values seem very out of scale with one another. But I don't know if they should stay that way, the article provided only the value of the spring constant. I am confused about that. I understood your remark about the constants, I agree, I had not thought about the rigidity due to the high values. Thank you so much for showing me this. I modeled the human operator (HO) as a mass+spring system. But for both HO and Haptic Interface, I couldn't define values for the constants, so I tried to follow the same logic.

• Welcome to Robotics, Leonardo Fagundes. The article you mention (and the figures you cite) are behind a pay wall. Can you post a copy of the figures here? That said, your damping and spring constant values seem very out of scale with one another. First I would never expect a negative damping (or else you'd get runaway acceleration) but also your spring constants are SO large as to be effectively rigid. Even without seeing your source figure, if you're guessing at terms, I would go with something like 10/1/2 for a spring/mass/damper system. But again, please show us what you're trying to do :) Apr 28 '20 at 2:51
• Hi, Chuck. Thank you very much :D I apologize, part of the paper can be viewed in ReseachGate by the link: link I posted a copy of the figures I mentioned, as you requested :) Apr 28 '20 at 18:18
• The published article follows in the link "Time-Domain Passivity Control of Haptic Interfaces", Hanaford and Ryu, 2002, IEEE transactions on robotics and automation vol .18, No1 :) Apr 28 '20 at 18:40
• robotics.stackexchange.com/questions/12197/… Apr 29 '20 at 22:11
• Thanks, @hauptmech But I've seen this post and I haven't been able to get the answer I'm looking for regarding the entry of the system mainly, how to insert the time delay... I also contacted the author, but got no response :( May 2 '20 at 17:48