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.
