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.
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:
Nonlinearity is time-invariant.
Nonlinearity does not generate any subharmonic as a response to the input sinusoidal.
The system filters out the super-harmonics generated by the nonlinearity (this assumption is often referred to as filtering hypothesis).
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: