Code covered by the BSD License  

Highlights from
Hydraulic Valve Parameters From Data Sheets and Experimental Data

image thumbnail

Hydraulic Valve Parameters From Data Sheets and Experimental Data

by

 

16 Apr 2010 (Updated )

Models and white paper on obtaining realistic parameter values from data sheets and measured data.

...
function F  = ...
    check_frequency_response(x,frequency_20,frequency_100, method)
% Checks the result of optimization by plotting phase characteristics at
% optimal parameters and computing the frequencies at -pi/2 phase shift.
% The frequency responses of a nonlinear system are obtained by processing
% the pulse transient characteristics with the FFT algorithm. 
% Copyright 2010 MathWorks, Inc.

% frequency_20   - frequency (Hz) at phase shift in -pi/2 at 20% input signal
% frequency_100  - frequency (Hz) at phase shift in -pi/2 at 100% input signal
% method - method to check the results, can be set to 'FFT' or 'direct'.

if strcmpi(method, 'FFT')
    model = 'actuator_freq_testrig_pulse_FFT_method';
elseif strcmpi(method, 'direct')
    model = 'actuator_freq_testrig_direct_method';
else
    error('Unspecified method');
end

load_system(model);

assignin('base','gain', x(1));
assignin('base','time_const', x(2));
assignin('base','saturation', x(3));

sim(model);

if strcmpi(method, 'FFT')
    
    y_20 = yout(:,2);               % Pulse transient characteristic at 20% input
    y_100 = yout(:,1);              % Pulse transient characteristic at 100% input
    fs = 1000;                      % Sampling frequency
    n = length(y_20);               % Window length = Transform length
    y_20_fft = fft(y_20,n);         % Discrete Fourier Transform
    y_100_fft = fft(y_100,n);       % Discrete Fourier Transform
    f0 = (0:n/2-1)*(fs/n);          % Shifted frequency range, positive region
    y_20_0 = fftshift(y_20_fft);    % Shifted DFT at 20% input
    y_100_0 = fftshift(y_100_fft);  % Shifted DFT at 100% input
    % Phase characteristic at 20% input for positive frequencies after unwrap
    phase_20 = unwrap(angle(y_20_0(257:end))); 
    % Phase characteristic at 100% input for positive frequencies after unwrap
    phase_100 = unwrap(angle(y_100_0(257:end)));

    % Computing frequency at 90 deg phase shift by interpolation of phase
    % characteristics
    frq_20 = interp1(phase_20,f0,-pi/2);
    frq_100 = interp1(phase_100,f0,-pi/2);


    % Computing frequency at phase shift in -pi/2 by interpolation of phase
    % characteristics
    frq_20 = interp1(phase_20,f0,-pi/2);
    frq_100 = interp1(phase_100,f0,-pi/2);
    % Computing errors
    error_20 = (frequency_20 - frq_20)/frequency_20 * 100;
    error_100 = (frequency_100 - frq_100)/frequency_100 * 100;


    disp(['******************************************************************']);  
    disp(['     System Characteristics with FFT at Obtained Parameters']);
    disp(['     Frequencies at -90 deg phase shift and 20% input:']);
    disp(['Target:  ', num2str(frequency_20),' Hz', ...
        '  Found:  ',num2str(frq_20),' Hz', ...
        '   Error :   ',num2str(error_20),'%']);
    disp(['     Frequencies at -90 deg phase shift and 100% input']);
    disp(['Target:  ', num2str(frequency_100),' Hz', ...
        ' Found:     ',num2str(frq_100),' Hz', ...
        '   Error:   ',num2str(error_100),'%']);

    
    % Plotting the phase frequency characterictics
    semilogx(f0,phase_20*180/pi,'LineWidth',3);
    hold on
    semilogx(f0,phase_100*180/pi,'r--','LineWidth',3), grid on;
    xlabel('Frequency (Hz)','FontSize',14)
    ylabel('Phase (degrees)','FontSize',14)
    title('Frequency Response (Phase)','FontWeight','Bold','FontSize',16)
    legend({'20% Input' '100% Input'},'FontSize',12);
    axis([2 40 -100 0]);
    plot(frq_100,-90,'ro','MarkerSize',10,'MarkerFaceColor','r')
    plot(frq_20,-90,'bo','MarkerSize',10,'MarkerFaceColor','b')
    hold off
    

else    % Direct measurement. yout contains input and output values of
        % both signals
        
    % Phase angle at specified frequency and input signal in 100%
    phase_100 = phase_computation(yout(:,[1,2]), tout, frequency_100,4);
    % Phase angle at specified frequency and input signal in 20%
    phase_20 = phase_computation(yout(:,[3,4]), tout, frequency_20,6);

    error_20 = (phase_20 + 90)/90 * 100;
    error_100 = (phase_100 +90)/90 * 100;


    disp(['******************************************************************']);  
    disp(['     System Characteristics after Direct Measurement at Obtained Parameters']);
    disp(['     Phase shift at ',num2str(frequency_20),'Hz',' and 20% input:']);
    disp(['Target:  ',' -90 deg', ...
        '  Found:  ',num2str(phase_20),' deg', ...
        '   Error :   ',num2str(error_20),'%']);
    disp(['     Phase shift at ',num2str(frequency_100),'Hz',' and 100% input']);
    disp(['Target:  ',' -90 deg', ...
        '  Found:   ',num2str(phase_100),' deg', ...
        '   Error:   ',num2str(error_100),'%']);
    disp(['******************************************************************']);  

end
end

function [phase] = phase_computation(yout, tout, freq_at_90, N)
% Function computes phase angle by processing sinusoids stored in yout

% yout(:,1)- input; yout(:,2) - output
% tout - simulation step time
% freq_at_90 - frequency at which phase angle equals -90 deg
% N - number of periods to settle the transients

% Measurement time. Time in N periods is assumed to be enough to settle
% down the transients. The input signal at this time is equal zero

t_s = 1/freq_at_90 * N;             % Start time
t_end = t_s + 1/freq_at_90;         % End time

% yout index at measurement time
for j = 1:length(tout)- 1
    if t_s >= tout(j) && t_s < tout(j+1)
        k_start = j;                % First array index after start time
        t_corr_1 = tout(j) - t_s;   % Correction due to discretization
    end
    if t_end >= tout(j) && t_end < tout(j+1)
        k_end = j-1;                  % Last array index within period
    end
end
% Computing time the output signal crosses zero
for j = k_start : k_end
    if yout(j,2) <= 0 && yout(j+1,2) > 0
        t_tab_cross = tout(j);   % First value in the table after crossing
        % Approximating real crossing time
        t_corr_2 = -yout(j,2) / (yout(j+1,2) - yout(j,2)) * ...
            (tout(j+1) - tout(j));
    end
end
% Crossing time for the output signal
t_cross = t_tab_cross + t_corr_1 + t_corr_2;
% Phase angle
phase = -(t_cross - t_s) / (1/freq_at_90) * 360;

end


% EOF

Contact us