Background on Describing Function Analysis

Describing function analysis is a widely known technique to study frequency response of nonlinear systems. It is an extension of linear frequency response analysis. In linear systems, transfer functions depend only on the frequency of the input signal. In nonlinear systems, when a specific class of input signal such as a sinusoidal is applied to a nonlinear element, you can represent the nonlinear element by a function that depend not only on frequency, but also on input amplitude. This function is referred to as a describing function. Describing function analysis has a wide area of applications from frequency response analysis to prediction of limit cycles.

To use sinusoidal-input describing function analysis, which is the most common type of describing function analysis, your model should satisfy these conditions:

In this example, you perform describing function analysis on a model with saturation nonlinearity that satisfies all the assumptions above.

Step 1: Opening the Model

Open the Simulink model with a saturation nonlinearity.

scdsaturationDF
mdl = 'scdsaturationDF';

Step 2: Describing Function Analysis of Saturation Nonlinearity

The saturation nonlinearity has the following sinusoidal input describing function:

where for a saturation with upper and lower limits of 0.5 and -0.5, respectively where A is the amplitude of the sinusoidal input signal.

Compute and plot the describing function, N_A(A), vs. amplitude, A, for amplitudes varying between 0.1 and 2.1:

A = linspace(0.1,2.1,21);
N_A = saturationDF(0.5./A);
plot(A, N_A);
xlabel('Amplitude');ylabel('N_A(A)');title('Describing function for saturation');

You can compute the describing function for saturation nonlinearity using FRESTIMATE over the same set of amplitudes for a fixed frequency of 5 rad/s. Note that the describing function of saturation does not depend on frequency, thus it is sufficient to run the analysis at a single frequency. Run a loop over all amplitudes where you will create a sinestream input with the (fixed) frequency and given amplitude, then run FRESTIMATE using this input signal at each iteration.

w = 5;
io(1) = linio('scdsaturationDF/In1',1,'input');
io(2) = linio('scdsaturationDF/Saturation',1,'output');
N_A_withfrest = zeros(size(N_A));
for ct = 1:numel(A)
    in = frest.Sinestream('Frequency',w,'Amplitude',A(ct));
    sysest = frestimate(mdl,in,io);
    N_A_withfrest(ct) = real(sysest.resp);
end
plot(A,N_A,A,N_A_withfrest,'r*')
xlabel('Amplitude');ylabel('N_A(A)');title('Describing function for saturation');
close_system(mdl);

Step 3: Closed-Loop Describing Function Analysis

You can also run a closed-loop describing function analysis over a frequency range. You can begin by analytically computing the frequency response from reference to output using describing functions. To do this, first compute the amplitude of the input signal "nonlinear_input" for the nonlinearity, given the reference amplitude and frequency. Note that the input amplitude for the nonlinearity is not necessarily equal to reference amplitude.

bdclose('scdsaturationDF');
scdsaturationDFcl
mdl = 'scdsaturationDFcl';

L = zpk([],[0 -1 -10],1);
w = logspace(-2,2,50);
A_DF = zeros(numel(A),numel(w));
for ct_amp = 1:numel(A)
    for ct_freq = 1:numel(w)
        % Compute the amplitude to nonlinearity solving the analytical
        % equation
        A_DF(ct_amp,ct_freq) = fzero(@(A_DF) solveForSatAmp(A_DF,L,w(ct_freq),A(ct_amp)),A(ct_amp),...
            optimset('Display','off'));
    end
end

Next, compute the analytical frequency response of the closed loop from reference to output with describing function for each amplitude and store it in an FRD-array.

L_w = freqresp(L,w);
for ct = 1:numel(A)
    N_A = saturationDF(0.5./A_DF(ct,:));
    cl_resp = N_A(:).*L_w(:)./(1+N_A(:).*L_w(:));
    cl(1,1,ct) = frd(cl_resp,w);
end

You can obtain frequency response for the closed loop from reference to input using FRESTIMATE in a similar way to the describing function analysis of saturation above.

io(1) = linio('scdsaturationDFcl/r',1,'input');
io(2) = linio('scdsaturationDFcl/Linear',1,'output');
for ct = 1:numel(A)
    in = frest.Sinestream('Frequency',w,'Amplitude',A(ct),...
        'NumPeriods',10,'SettlingPeriods',7);
    cl_withfrest(1,1,ct)  = frestimate(mdl,in,io);
end

You can plot analytically calculated closed-loop magnitude along with the one from FRESTIMATE.

h = figure;bodemag(cl,'b',cl_withfrest,'r');
annotation(h,'textarrow',[0.64 0.58],[0.64 0.58],'String','Increasing A');

Close the model:

bdclose('scdsaturationDFcl')