Describing Function Analysis of Nonlinear Simulink Models
This example shows how to use the frequency response estimation to perform a sinusoidal-input describing function analysis for a model with a saturation nonlinearity.
Describing function analysis is a technique for studying the 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 sinusoid, is applied to a nonlinear element, you can represent the nonlinear element using a describing function. A describing function depends not only on the input frequency but also on the input amplitude. Describing function analysis has a wide area of applications from frequency response analysis to the 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 the following conditions.
Nonlinearity is time-invariant.
Nonlinearity does not generate any subharmonics in response to the input sinusoid.
The system filters out any superharmonics generated by the nonlinearity (this condition is often referred to as the filtering hypothesis).
In this example, you perform describing function analysis on a model with a saturation nonlinearity that satisfies these conditions.
Open the Simulink® model.
mdl = 'scdsaturationDF'; open_system(mdl)
The saturation nonlinearity has the following sinusoidal input describing function.
Here, for a saturation with upper and lower limits of 0.5 and -0.5, respectively, and A is the amplitude of the sinusoidal input signal.
Compute and plot the describing function for amplitudes varying between 0.1 and 2.1. The describing function is implemented in
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 the saturation nonlinearity using
frestimate over the same set of amplitudes for a fixed frequency of 5 rad/s. Since the describing function for a saturation does not depend on frequency, it is sufficient to run the analysis at a single frequency.
Define the input and output linear analysis points.
io(1) = linio('scdsaturationDF/In1',1,'input'); io(2) = linio('scdsaturationDF/Saturation',1,'output');
For each amplitude, create a sinestream input with the fixed 5 rad/s frequency and given amplitude. Then, run
frestimate using this input signal.
w = 5; 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 the frequency response estimation result along with the describing function.
plot(A,N_A,A,N_A_withfrest,'r*') xlabel('Amplitude') ylabel('N_A(A)') title('Describing Function for Saturation Using frestimate') legend('Describing function','Frequency response estimate')
Close the model.
You can also run a closed-loop describing function analysis over a frequency range. Open the model for this analysis.
mdl = 'scdsaturationDFcl'; open_system(mdl)
To begin, analytically compute the frequency response from the reference to the output using describing functions. To do so, first compute the amplitude of the input signal to the nonlinearity, given the reference amplitude and frequency. The input amplitude for the nonlinearity is not necessarily equal to the reference amplitude.
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 input amplitude to the nonlinearity by solving the % analytical equation. opt = optimset('Display','off'); A_DF(ct_amp,ct_freq) = ... fzero(@(A_DF) solveForSatAmp(A_DF,L,w(ct_freq),A(ct_amp)),... A(ct_amp),opt); end end
Next, for each amplitude, compute the analytical frequency response of the closed-loop system from the reference to the output using the describing function. Store the results in an array of
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
Obtain the frequency response for the closed-loop system using
frestimate in a manner similar to the saturation describing function analysis.
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
Plot the analytically calculated closed-loop magnitude along with the result from
frestimate. As expected, the results match. The arrow indicates the direction of increasing amplitude.
h = figure; bodemag(cl,'b',cl_withfrest,'r') annotation(h,'textarrow',[0.64 0.58],[0.64 0.58],'String','Increasing A'); legend('Analytical result','Frequency response estimate',... 'Location','southwest')
Close the model.