Describing Function Recording with Simulink and MATLAB

Describing function is an equivalent gain of nonlinear element, defined by the harmonic linearization method of nonlinear static characteristic (Novogranov, 1986, Slotine and Li, 1991, Schwarz and Gran, 2001, Vukic et al. 2003 and many others). It is a known method of analysis and synthesis when nonlinear system can be decoupled into linear and nonlinear parts (Fig. 1). If the linear part of the system has the characteristics of low-pass filter and if we apply periodical signal to the system, output signal will have the same base frequency as input signal with damped higher frequencies.


Introduction
Describing function is an equivalent gain of nonlinear element, defined by the harmonic linearization method of nonlinear static characteristic (Novogranov, 1986, Slotine and Li, 1991, Schwarz and Gran, 2001, Vukic et al. 2003 and many others).It is a known method of analysis and synthesis when nonlinear system can be decoupled into linear and nonlinear parts (Fig. 1).If the linear part of the system has the characteristics of low-pass filter and if we apply periodical signal to the system, output signal will have the same base frequency as input signal with damped higher frequencies.If the amplitudes of higher harmonics are relatively small when compared to the amplitude of the first harmonic, i.e. if filter hypothesis holds, output signal can be approximated by it's first harmonic.Mathematically, first harmonic of the output signal, for the sinusoidal input signal sin m Xt  , can be expressed by the Fourier expressions: Describing function is the ratio between first harmonic of the output signal and input signal in complex form: where   m PX and   m QX are coefficients of the harmonic linearization (Novogranov, 1986, Slotine and Li, 1991, Schwarz and Gran, 2001, Vukic et al., 2003 and many others).
Determination of describing function boils down to the determination of integral expressions for the known static characteristic of the nonlinear part of the system.
There are many conventional nonlinearities for which static characteristics and describing functions are theoretically derived and given in analytical form.
The problem arises when the static characteristic of nonlinear system cannot be analytically expressed or integral expressions can not be solved.In that case describing function can be determined by experiment or simulation (Kuljaca et al., Sijak et al.k 2004, 2005, 2007a, 2007b) or some method of numerical integration (Schwarz and Gran, 2001).

Nonlinear elements in Simulink
Nonlinear elements are given in Simulink library Discontinuities.Twelve discontinuities given there are shown in Fig. 2. We can see that there are three elements called "Dynamic" (Dead Zone Dynamic, Saturation Dynamic and Rate Limiter Dynamic).The term dynamic is considered at these elements with respect to change of nonlinearity limits (output values limits or rate change limits), not as dynamics in control systems sense (i.e.behavior in time domain).Never the less, we will not deal with such elements here.Describing function method analysis or synthesis is not suitable for systems with varying parameters.
One can see that Simulink gives only basic non-dynamic nonlinearities.More nonlinearities and their describing function can be found in work by Vukic et al., (2003).It is clear that one needs to build them out of basic Simulink models or write them as m or s functions.Since given functions don't have dynamics it is recommended to build them out of given Simulink blocks, or if that is not possible, then using type 2 s-functions.Even in that case some real time extension software for Simulink and Matlab might not be able to handle user designed s-functions.

Describing function by numerical integration
An interesting approach to obtaining describing function using Simulink is given by Schwartz and Gran.The approach is very simple, yet effective.In essence, equations (2) and (3) are directly applied in Simulink on nonlinear element outputs.A simulation scheme adapted from Schwartz is given in Fig. 3.In this case dead zone nonlinearity is analyzed.Instead of dead zone, any other nonlinearity can be included in model.5).It can be seen that imaginary part of the gain is zero.That is a case with all symmetrical nonlinearities. - where A m is amplitude of signal entering saturation.
Let us now see the results of running our simulation script in order to obtain describing function for saturation.Simulation results are shown in Fig. 8 and Fig. 9.It can be seen from Fig. 8 that describing function in this case is not dependent on frequency.Also, phase is zero since saturation is symmetrical nonlinearity.That can be seen from very small values of describing function phase as given in Fig. 9.These values are practically zero.Obtained result from simulation is given complex matrix s.In case like this, where there is no frequency dependency and phase is zero, there is no need for plotting describing function in three dimensions.To make function plot in two dimensions we run the following short script: We will see in subsequent sections that this representation of symmetrical nonlinearities as gain dependent on amplitude of nonlinearity input signal is important for stability analysis.
This technique can be used to obtain describing functions of real nonlinearities by using Simulink model from Fig. 3 with appropriate real time additions with Real Time Window Target, xPc Target or other real time systems for use with Simulink.Of course, in that case amplitudes will have to be given one by one, not as here as one vector.

Fuzzy control systems and their describing functions in Simulink
Use of fuzzy control systems is today well known in control systems design.Unfortunately, while there are well developed methods for design or training of fuzzy systems, stability analysis remains problematic.If the fuzzy system is of Mamdani type, than its mathematical description can be very complicated and stability analysis is complicated as well.However, if we can obtain describing function of fuzzy system, then we can use known stability analysis methods developed for linear systems.One type of experimental method of obtaining describing function for fuzzy elements is given in (Kuljaca et al.) and (Sijak et al., 2005(Sijak et al., , 2007b)).
Here, we will use numerical integration in Simulink as given by Schwartz.Simulation file will have to be reworked since we cannot use vectors in this simulation due to fuzzy element.
Of course, there are different fuzzy controllers and elements, but describing function for all of them can be achieved using this method.One should note that here we are not talking about adaptive fuzzy controllers, but about fuzzy controllers with fixed membership functions and rule base.
Let us first look on one example of fuzzy controller, as shown in Fig. 11   First of all, it can be seen that fuzzy element here has two inputs.Fuzzy elements in general can have more inputs and outputs than here, but, this structure is most often in control systems (either as given here with derivative of the error or with integral of the error).When more inputs are used it becomes extremely complicated to tune the controller.
Numerical integration Simulink model (Schwartz) with fuzzy controller will look like in Fig. 12, without vector inputs for amplitudes since fuzzy block in Simulink cannot handle such input type.This is also much closer to real measurement, since in real measurement we would not able to use vector inputs.Meaning of this is that we need to use adjusted script to run the Simulink model.Adjusted script is given in Fig. 13.
Fuzzy Logic controller Simulink block is regular block from Simulink Fuzzy Logic toolbox library.
Fuzzy element itself is developed using FIS editor from Matlab Fuzzy Logic Toolbox.This is an illustrative example; however, it is quite useful in giving an insight in use of describing function for fuzzy controllers.The given method can give a graphic representation of any fuzzy controller based on error and derivatives of error signal.In most cases fuzzy controllers are also symmetrical, thus subsequent stability analysis will be simpler.Amplitudes and frequencies are to be chosen to satisfy expected operational environment of controller.In example given here gains k p and k d are set to 1.2 and 0.001 respectively, min-max and centroid defuzzification Mamdani type fuzzy controller membership functions are given in Fig. 14 and rulebase is given in Fig. 15.Describing function phase as function of amplitude and frequency shown in Fig. 17 shows very small values and for all practical purposes can be considered to be zero.Since our fuzzy controller design is symmetrical one that is in compliance with theoretical expectations.

Stability analysis with describing functions in Simulink -Case study
Once describing function is obtained, it could be used for stability analysis.We will use a simplified secondary frequency control model of the isolated thermo power system with generation rate constraint (Sijak et al.,2001).Control model is shown in Fig. 18. -the steam turbine time constant P m -the change of the mechanical power of the turbine P L -the change of the power system load P r -the power system active power reference change f -the power system frequency change -the power system transfer function R -the static speed droop of the uncontrolled system () Fz -the power system generation rate constraint, static characteristic of saturation nonlinearity.
The system given in Fig. 18 consists of nonlinear parts divided by linear parts for which the filter hypothesis is satisfied.Such system can be represented as in Fig. 19.

 
Li G p , the nonlinearities can be harmonically linearized and their describing functions used instead (Vukic et al., Netushil et al.).
Harmonically linearized system is shown in Fig. 20. and nonlinear part of the plant respectively, then the stability analysis of the system can be conducted using the solution of complex equation (7).
Let us now deal with nonlinearities.Generation rate constraint is saturation type nonlinearity shown in Fig. 21.Describing function for such nonlinearities is known and it is given in (8).Fuzzy controller represents a bit more complex problem.We do not have its describing function in analytical form and we need to obtain it by experimental method.The method is described in Section 5 and describing function gain is plotted in Fig. 16 in three dimensional representation.Phase changes are zero, thus there is no complex component of describing function.
Representation in three dimension in this case is not required since there is no dependency of describing function gain on frequency.After running simulation script given in Section 5, one should run the following Matlab code in order to obtain two-dimensional representation of fuzzy controller describing function:  The stability of the equilibrium point for the system given in Fig. 18 can be determined from the characteristic polynomial of the closed loop harmonically linearised equation (Kuljaca et al., Sijak et al., 2007): where: 0 1 K R  From (9) the closed loop characteristic equation of the controlled system can be obtained: The system is stable as long as ( 11) is positive.Thus, the stability boundary can be derived as:  With given values system is stable in operating region of interest.

Conclusion
Chapter deals with use of Simulink for obtaining describing function for different elements.Described method of numerical simulation performed in Simulink is suitable for obtaining describing function of any single input single output non-dynamic nonlinearity.Method used for conventional static nonlinearities is then adjusted to work with fuzzy systems (and other systems that cannot handle vector based batch simulation in Simulink).Matlab scripts required to run given Simulink model are also given in Chapter.This method can be also used to obtain describing function for single input single output static neural networks.Finally, an example from use describing function for stability analysis is given.
Described method can be used with practical controllers with Matlab real time tools.

Fig. 3 .
Fig. 3. Simulink model 'csm' for generating describing function (Schwartz) Model is called from Matlab script (adapted form Schwartz and Gran as well) as shown in Fig. 4. simulation parameters are passed to Simulink schematic from Matlab script, as given in Fig. 5.

Fig. 4 .
Fig. 4. Matlab script for describing function generation Fig. 8. Describing function gain as function of amplitude and frequency

Fig. 9 .
Fig. 9. Describing function phase as function of amplitude and frequency Fig. 10.Describing function gain -two dimensional representation (Kuljaca et al., Sijak  et al. 2007a), where k p and k d are gains.

Fig. 14 .
Fig. 14.Membership functions of the fuzzy element: a) proportional input, b) derivative input and c) output . It can be seen from Fig.16that describing function gain is function of amplitude and it really does not depend on frequency in this case.
Fig. 16.Describing function gain as function of amplitude and frequency

Fig. 17 .
Fig. 17.Describing function phase as function of amplitude and frequency Fig. 18.Power system secondary load-frequency control model where: 1 1 G G G sT   -the transfer function of the turbine governor CH T-the steam turbine time constant P m -the change of the mechanical power of the turbine P L -the change of the power system load P r -the power system active power reference change f -the power system frequency change

Fig. 19 .
Fig. 19.The structure of the nonlinear system and fuzzy controller where:   ,, ii i d Fxp x p dt  , -nonlinear parts of the system  Li G s -transfer functions of the linear parts of the system Assuming that the filter hypothesis is valid for

Fig. 20 .FX
Fig. 20.The structure of the harmonically linearised nonlinear system and fuzzy controller The self-oscillations (limit cycles) of the system in Fig. 20 are described by solution of:        11 1 22 2 1, , 0 mL m L FX G j FX G j    (7)Assuming that   11 , m FX  and   22 , m FX  are describing functions of the fuzzy controller Fig. 21.Generation rate constraint gdf=s(:,1); plot(E,real(gdf)),grid, ylabel('Describing function gain'); xlabel('Amplitude');Describing function in two dimensional representation is given in Fig.23.

Fig. 23 .
Fig. 22.Generation rate constraint describing function 10)By applying Hurwitz stability criterion the following inequality is obtained: used for simulation, the stability boundary function G NF = f(G NZ ) can be evaluated.The function G NF = f(G NZ ) is shown in Fig. 24.