Open access peer-reviewed compact

# Rational Fitting Techniques for the Modeling of Electric Power Components and Systems Using MATLAB Environment

By Eduardo Salvador Bañuelos-Cabral, José Alberto Gutiérrez-Robles and Bjørn Gustavsen

Submitted: May 26th 2017Reviewed: September 29th 2017Published: December 6th 2017

DOI: 10.5772/intechopen.71358

## Abstract

This book provides a detailed description of some of the most widely used rational fitting techniques for approximation of frequency domain responses. The techniques are: Bode’s asymptotic approximation, the Levy method, iteratively reweighted least squares, the Sanathanan-Koerner method, the Noda method, Vector Fitting, the Levenberg-Marquardt method, and the Damped Gauss-Newton method. A MATLAB routine for each technique is presented. These techniques are tested by approximating synthetic frequency domain responses. Then, they are applied to the rational approximation of the frequency-dependent parameters corresponding to a single-phase transmission line. The effect of the rational function-based models is evaluated, considering transients in three cases: Open-ended, short-circuited, and perfectly matched lines. The error levels obtained in time domain simulations are consistent with the fitting deviations of the frequency-dependent parameters. The book concludes by showing main advantages and disadvantages for each technique.

### Keywords

• rational approximation
• least squares
• weighted least squares
• Bode’s asymptotic approximation
• vector fitting

## 1. Introduction

One of the main problems encountered in the modeling of power system components is the inclusion of frequency-dependent effects in a time domain simulation . In practice, these effects are described by discrete frequency domain responses that are obtained via calculations or measurements . For the time domain simulation, implementation of numerical convolutions is possible but is computationally inefficient. The frequency-dependent effects are usually performed in the frequency domain via complex-curve fitting processes, leading to rational function-based models which can be expressed in pole-zero form, pole-residue form, and polynomial form . Once obtained, these models must be converted into a form suitable for representation in an electromagnetic Transient Program (EMTP-type) circuit simulator . Commonly, such representations can be in the form of network synthesis or in ordinary differential equations (ODEs). Some methods have been adapted to particular problems due to the difficulty present in the development of a general methodology. These rational fitting techniques have been developed since the 1950s for model synthesis based on frequency response data. Some antecedents are presented below.

In 1959, Levy  presented a mathematical procedure of linearized Least Squares (LS) for model synthesis based on a polynomial form. Some years later, in 1963, Sanathanan and Koerner  improved the Levy method by introducing an iterative weighted method to reduce the biasing in the approximation caused by the linearization. Also, based on the Levy method, Lawrence and Rogers  presented in 1979 a sequential algorithm that allows the fitting progress of a transfer function point-to-point without matrix inversion.

In addition, a technique based on pole-zero form that quickly became popular for the modeling of overhead lines was presented by Marti  in 1982. This method is a numerical implementation of the well-known Bode diagrams technique. Recently, in 2017, Marti and Tavighi  presented an investigation based in the same rational approximation technique for the modeling of transmission lines. In this work, this fitting technique will be referred to as Asymptotic Approximation or Bode.

In 1993, Soysal and Semlyen  used the Gauss–Seidel optimization method to improve the results given by Levy and, 6 years later Gustavsen and Semlyen  developed the Vector Fitting method which has become one of the most widely used techniques.

In 2006, Gustavsen  proposed a modification of his algorithm in order to improve the ability of VF to relocate poles to better positions. This is achieved by replacing the high-frequency asymptotic requirement of the VF scaling function with a more relaxed condition. He called this method Relaxed VF (RVF).

The VF method is closely related to the Universal Line Model (ULM) . The ULM is formulated in terms of rational approximation of the line parameters through VF. A recent publication, in 2016,  Bañuelos et al. propose the use of only real poles and zeroes in the ULM to improve its numerical efficiency.

In 2005, Noda  presented an iterative algorithm that partitions the entire frequency range in order to avoid ill-conditioning of the system when Levy method is used. Recently, in 2015, Bañuelos-Cabral et al.  propose the implementation of Damped Gauss-Newton (DGN) to increase the accuracy of this technique.

The aim of this book is to provide a MATLAB algorithm and a detailed description of the rational approximation techniques that are most commonly used to approximate functions in frequency domain. These techniques are: Bode’s Asymptotic Approximation (Bode), the Levy (Levy or LS) method, Iteratively Reweighted Least Squares (IRLS), the Sanathanan-Koerner (SK) method, the Noda (Noda) method, Vector Fitting (VF), the Levenberg–Marquardt (LM) method and the Damped Gauss-Newton (DGN) method.

In this chapter, the methodology for determining the state-space (SS) representation in concordance with the model used in the approximation is presented first. Next, the abovementioned fitting techniques are described in detail. Then, the techniques are tested by approximating a synthetic frequency domain response. The techniques are implemented in MATLAB environment. Finally, advantages of the methods are demonstrated in the rational approximation of the frequency-dependent parameters corresponding to a single-phase transmission line.

## 2. State-space model from a transfer function

The SS representation of a given system can be determined from its frequency response on rational form. Basically, there are three different rational forms (models) that can represent a measured or calculated frequency response: pole-zero form, polynomial form, and pole-residue form. In this section, a methodology is presented for determining the SS representation according to the given model type.

### 2.1. Pole-zero form

Usage of pole-zero form or series realization leads to a rational function-based model of nth-order given by (1), which is a ratio between products of first-order transfer functions,

Fsksz1sz2sznsp1sp2spm.E1

The generalized graphic representation of pole-zero form with n = m is shown in Figure 1 .

From Figure 1, it is possible to obtain the state equations and the output equation as

ẋ1=p1x1+uẋ2=p2x2+z1x1+ẋ1ẋ3=p3x3+z2x2+ẋ2ẋn1=pn1xn1+zn2xn2+ẋn2ẋn=pnxn+zn1xn1+ẋn1E2
y=kznxn+ẋn.E3

Using algebraic manipulation in (2) and (3) to remove ẋnfrom the right side, the state equations and output equation in matrix form are

ẋ1ẋ2ẋn1ẋn=p1000ĉ1p200ĉ1ĉ2pn10ĉ1ĉ2ĉn1pnx1x2xn1xn+1111uE4
y=ĉ1ĉ2ĉn1ĉnx1x2xn1xn+kuE5

where ĉn=kznpn. There is a direct relation between the SS coefficients in (4) and (5) with the transfer function coefficients in (1) for pole-zero form realization. Thus, the SS description can be obtained directly from the transfer function by inspection.

### 2.2. Polynomial form

In this case, polynomial form or direct form realization leads to a rational function-based model of nth-order given by (6), which is a ratio of two polynomials,

Fsa0+a1s+a2s2++ansnb0+b1s+b2s2++bmsm.E6

A graphical representation of the polynomial form with bm = 1 and n = m is shown in Figure 2 .

From Figure 2, it is possible to obtain the state equations and output equation as

ẋ1=x2ẋ2=x3ẋn1=xnẋn=b0x1b1x2bn1xn+uE7
y=a0x1+a1x2++an1xn+anẋn.E8

To remove ẋnfrom (8) the last equation of (7) can be used. These equations can be written more conveniently as

ẋ1ẋ2ẋn1ẋn=010000010000001b0b1b2bn2bn1x1x2xn1xn+0001uE9
y=â0â1ân1x1x2xn+anu,E10

where âi=aianbi. Also, there is a direct relation between the SS coefficients in (9) and (10) with the transfer function coefficients in (6). The matrices of the SS description can be obtained directly from the transfer function through inspection.

### 2.3. Pole-residue form

Finally, pole-residue form or parallel realization leads to rational function-based model of nth-order given by (11), which is a sum of partial fractions

Fsc1sp1+c2sp2++cnspn+d.E11

The generalized graphic representation of the pole-residue form is shown in Figure 3 .

From this Figure 3, it is possible to obtain the state equations and output equation as

ẋ1=p1x1+uẋ2=p2x2+uẋn=pnxn+uE12
y=c1x1+c2x2++cnxn+du.E13

In matrix form it yields

ẋ1ẋ2ẋn=p1000p2000pnx1x2xn+111uE14
y=c1c2cnx1x2xn+du.E15

Similarly to the pole-zero form and the polynomial form, the SS representation can be obtained through inspection. It should be mentioned that the pole-residue form produces an uncoupled SS system.

## 3. Fitting methods

This section provides a description of the most widely used techniques for rational approximation of frequency domain responses within the power systems area.

### 3.1. Asymptotic approximation (Bode) technique

The Asymptotic Approximation technique consists of a numerical implementation of the graphical technique known as Bode diagrams . The Asymptotic Approximation was first implemented by Marti  to include the frequency dependence in transmission line modeling. This technique considers the pole-zero form with real poles and zeros only,

Fs=ksz1sz2sznsp1sp2spm.E16

Considering s = , we get

Fs=kz1z2znp1p2pm.E17

Eq. (17) can be expressed in a standard form as

F=kz11z21zn1z1z2znp11p21pm1p1p2pmE18

and in polar form

F=K0z11α1z21α2zn1αnp11β1p21β2pm1βmE19

where K0 = kz1z2zn/p1p2pm. The Asymptotic Approximation technique approximates only the magnitude of the frequency response, which can be expressed as

F=K0z11z21zn1p11p21pm1.E20

Finally, by using logarithmic properties in Eq. (20), it is obtained:

log10F=log10K0+log10/z11++log10/zn1log10/p11log10/pm1.E21

In Bode diagrams, each term of (21) is plotted individually in accordance with its already known asymptotic behavior and combined in order to obtain the desired diagram.

In the numerical implementation, the function to be fitted is compared with the sum of the line segments given by the addition of a pole or a zero in the model (21); the precision of the fitting depends on the sensitivity to locate these poles and zeros into the model.

### 3.2. Levy (Levy) method

This method was introduced by Levy  for complex-curve fitting. The Levy method makes the identification of the polynomial coefficients (22) in a LS sense. It is also known simply as Lest Squares (LS).

FsNsDs=a0+a1s+a2s2++ansn1+b1s+b2s2++bmsm.E22

The numerical difference between the frequency response to be fitted and the model represents the error in the approximation, that is,

εs=FsNsDs.E23

Multiplying both sides of Eq. (23) by D(s), we obtain

εs=εsDs=FsDsNs.E24

Considering that ε′ tends to zero, (24) can be expressed as

Fs1+b1s+b2s2++bmsma0+a1s+a2s2++ansn=0.E25

The unknown coefficients in (22) can now be solved by formulating (25) as a LS problem (Ax = b). Frequently, the sample size of the frequency response is greater than the number of polynomial coefficients to be calculated, so an overdetermined system is obtained, as follows:

Ak=1sksknskFskskmFskE26
x=a0a1anb1b2bmTE27
b=Fs1FskT,E28

where k denotes the k-th data sample. The objective function to be minimized is

minxAxb2.E29

This is called a linear least squares problem. The solution x satisfies the normal equations :

ATAx=ATb.E30

Finally, the LS solution is obtained by

x=ATA1ATb.E31

### 3.3. Weighted least squares (WLS)

The LS solution of a given system equation assumes that each equation (row in A) is equally important. However, the system equation itself may be biased, e. g. in the Levy method which is biased due to the multiplication of F(s) with D(s) in (24). In the rational approximation of frequency domain responses, Weighted Least Squares (WLS) is used to mitigate the biasing by giving more weight to specific equations in order to overcome the own deficiency of the used technique. Iteratively Reweighted Least Squares (IRLS), the Sanathanan-Koerner (SK) method and the Noda (Noda) method are techniques that implement the concept of WLS. These techniques are described below.

#### 3.3.1. Iteratively reweighted least squares (IRLS) iteration

A robust regression procedure is an alternative to LS solution when data are contaminated with outliers or influential observations (measurement and/or computation errors) [16, 17]. The idea of robust regression is to weight (less) these observations differently based on the proposal of weighting functions. IRLS can be considered as a robust regression procedure.

It is proposed to use IRLS for the rational approximation of frequency domain responses by using the weighting functions in inverse form, namely, data are not considered to be contaminated with outliers or influential observations; the error in the fitting in each iteration is used to weight (more) the frequency response data that have not been approximated correctly.

In WLS, the objective function to be minimized is

minxWbAx2E32

where W is a diagonal weighting matrix. Eq. (32) is called a linear weighted least squares problem and x the linear weighted least squares solution of the system. This solution satisfies the normal equations :

ATWAx=ATWb.E33

Then, the WLS solution is given by

x=ATWA1ATWb.E34

It is possible to solve (34) through an iterative process (35), where i is the iteration number.

xi+1=ATWiA1ATWib.E35

Eq. (35) is the IRLS method. According to (35), W must be updated iteratively; the error in the approximation is used for this purpose, which can be expressed as

εsx=FsF̂sx,E36

where F̂sxis the model to be fitted. Each element of W is updated according to

wki=wki1ψεki1,E37

with ψ(ε) being the weighting function. In the beginning of the process, W is an identity matrix.

Table 1 shows a list of weighting functions proposed to use in the implementation of IRLS. Weighting functions 1, 2, and 3 are the inverse of the functions used in robust regression [16, 17], and functions 4, 5, and 6 are proposed in this work.

Numberψ(ε)
1|ε|
21+ε21
31εp2,p=1.2
4ε21+ε2
5ε21+ε2
6εpp,p=1.2

### Table 1.

Weighting functions.

Ultimately, in Figure 4 the behavior of the weighting functions with ε = [−2, 2] is shown.

#### 3.3.2. Sanathanan-Koerner (SK) iteration

This method is based on the polynomial form:

FsNsDs=a0+a1s+a2s2++ansn1+b1s+b2s2++bmsm.E38

The objective is to identify the coefficients for the polynomials N(s) and D(s) in (38) that minimize the error function

εs=FsDsNs.E39

Sanathanan and Koerner  proposed an iterative procedure where (39) is divided by the denominator from the previous iteration; thus, (39) can be expressed as

εs=1Dsl1FsDslNsl.E40

The subscript l denotes the iteration number, and D(s) is considered 1 in the first iteration.

#### 3.3.3. Noda (Noda) iteration

This algorithm proposed by Noda  for the rational fitting identification is based on the function:

FsNsDs=a0+a1s+a2s2++ansn1+b1s+b2s2++bmsm.E41

This technique corresponds to IRLS using the weighting function number 1 in Table 1. Additionally, Noda also proposed a procedure to prevent the ill-conditioning of the system by partitioning the given frequency response data into sections along the frequency axis. The technique is then applied to each section of the frequency response in order to identify the poles, and finally, the corresponding residues are obtained by means of a standard LS procedure, using the entire frequency response.

### 3.4. Vector Fitting (VF) method

Vector Fitting performs the fitting by replacing a set of heuristically calculated initial poles with a set of relocated ones through an iterative procedure . VF works in two stages: First, it improves the position of the initial poles iteratively. Second, it calculates the residues in one step. This method is based on the pole-zero form

Fsn=1Ncnspn+d+sh,E42

where cn are the residues, pn are the poles, d is the constant term and h the proportional part.

#### 3.4.1. Pole identification

Considering the initial poles as an, and multiplying F(s) by an auxiliary function σ(s) gives

σsFsσFfits=n=1Ncnsan+d+shE43

where σ(s) is defined as

σsσfits=n=1Nc˜nsan+1,E44

with σFfit(s) being the fitting of σ(s)F(s) and σfit(s) being the fitting of σ(s). From (43) one can obtain,

FsσFfitsσs.E45

Substituting (43) and (44) into (45) we obtain,

Fshn=1N+1szn/n=1Nsann=1Nsz˜n/n=1Nsan=hn=1N+1sznn=1Nsz˜n.E46

Eq. (46) indicates that the zeros of σ(s) are an approximation of the poles of F(s). To obtain the zeros one multiplies (44) by F(s) which results in

σsFsn=1Nc˜nsan+1Fs.E47

Equating Eqs. (43) and (47) yields

n=1Ncnsan+d+sh=n=1Nc˜nsan+1Fs.E48

Algebraic manipulation in (48) gives

n=1Ncnsan+d+shn=1Nc˜nFssan=Fs.E49

Eq. (49) can now be formulated as a LS problem (Ax = b). Usually, the sample size of the frequency response is greater than the number of coefficients to be calculated, so an overdetermined system is obtained:

Ak=1ska11skaN1skFskska1FskskaNE50
x=c1cNdhc˜1c˜NTE51
b=Fs1FskTE52

The LS solution of the system (Ax = b) delivers the residues of σfit(s); the zeros of this function are the poles an for the next iteration, which correspond to the eigenvalues of ,

H=Gkc˜T,E53

where G is a diagonal matrix containing the poles and k is a unit column vector.

#### 3.4.2. Residue identification

Once the poles pn for the rational approximation in (42) have been obtained, the residues can be found by solving (42) as a LS problem. The new overdetermined system is solved similarly as (50), (51) and (52).

#### 3.4.3. Relaxation (RVF)

The scaling function σ(s) (44) is observed to approach unity at high frequencies. This asymmetry with respect to the frequency gives a tendency to relocate poles in the direction of lower frequencies. In addition to this biasing, it also reduces the convergence speed and often leads to a reduced accuracy for the final model. This problem is effectively alleviated by the introduction of Relaxed Vector Fitting (RVF) , where the scaling function (44) is replaced by

σsσfits=n=1Nc˜nsan+d˜E54

with d˜as a real variable. To obtain a non-trivial solution, a single line is added to the LS problem, where the sum of the real part of σ(s) over the frequency samples is fixed. This relaxed criterion allows σ(s) to freely vary in shape.

### 3.5. Levenberg-Marquardt (LM) method

The Levenberg–Marquardt (LM) method is a technique used to improve iteratively parameter values in order to reduce the sum of the squares of the errors between the model and the calculated or measured data [18, 19]. This method is actually a combination of two minimization methods: The gradient descent (GD) method and the Gauss-Newton (GN) method.

Consider the model F̂sx, which is a function of s and n unknown parameters contained in vector x. The model is to be fitted to a set of frequency response data F(s) as follows:

FsF̂sx,E55

where x = [x1, x2, ⋯xn] is a vector that contains the coefficients to be fitted. The LS error between the fitting of the model and the frequency response data is

γx=12k=1NsF̂skxFsk2E56
γx=12F̂sxFsTF̂sxFsE57
γx=12F̂sxTF̂sxF̂sxTFs+12FsTFs.E58

The aim of LM method is to find iteratively a perturbation h to modify the parameters x that reduces the objective function γ(x) .

#### 3.5.1. Gradient descent (GD) method

The gradient of a function indicates the direction of the maximum rate of change in a specific point of the function. GD updates the parameter values in the opposite direction to the gradient of the objective function. Thus, the gradient of γ(x) is:

xγx=F̂sxFsTxF̂sxFsE59
xγx=F̂sxFsTJ,E60

where J is the Jacobian matrix which represents the local sensitivity of the model F̂sxto variation in the parameters x. Finally, the perturbation h that moves the parameters in the direction of steepest descent is given by

hGD=JTF̂sxFsE61

#### 3.5.2. Gauss-Newton (GN) method

The Gauss-Newton (GN) method supposes that the objective function γ(x) is approximately quadratic in the parameters near the optimal solution. Approximating the model with a first-order Taylor series,

F̂sx+hF̂sx+F̂sxxh=F̂sx+Jh.E62

Substituting (62) into (57) gives

γx+h=12F̂TF̂F̂TF+12FTF+F̂FTJh+12hTJTJh.E63

The perturbation h that minimizes γ(x) is reached when

γx+hh=0.E64

Taking the derivative for (63) gives

γx+hh=F̂sxFsTJ+hTJTJE65

and applying condition (64), we obtain the perturbation h calculated by the GN method

hGN=JTJ1JTF̂sxFs.E66

#### 3.5.3. Levenberg-Marquardt method (LM)

LM works more like the GD method when the parameters are far from their optimal value, and more like the GN method when the parameters are close to their optimal value. Levenberg  and Marquardt  proposed this method:

JTJ+λIhLM=JTF̂sxFs,E67

where λ is known as the Levenberg–Marquardt parameter. Clearly, for large values of λ, the LM method results in a GD update; large distance from the function. The GD method is utilized to provide steady and convergent progress toward the solution. As the solution approaches the minimum, λ is adaptively decreased; then, the LM method approaches the GN method, and the solution typically converges rapidly to the local minimum.

The following modification has also been suggested:

JTJ+λdiagJTJhLM=JTF̂sxFsE68

There are different methods for update the Levenberg-Marquardt parameter λ. A complete review is presented in Ref. . Nevertheless, a simple method consists of reducing this parameter as the solution approaches the optimal value. The numerical implementation consists iteratively of:

1. Calculate the objective function γ(x), (57).

2. Calculate the Jacobian J, taking into account the implemented model (1), (6) or (11).

3. Calculate the perturbation h, (67).

4. If γ(x + h) ≤ γ(x), then x is replaced by x + h and λ is reduced.

5. If γ(x + h) > γ(x), then x is not replaced and λ is increased.

6. Convergence is achieved when the parameters reach a tolerance value or when the iteration count exceeds a pre-specified limit.

#### 3.5.4. Damped gauss-Newton (DGN)

The DGN method can be applied to any of the alternative function-based models: pole-zero form (1), polynomial form (6), and pole-residue form (11). In the DGN method, a damping factor α is introduced

hDGN=αJTJ1JTF̂sxFs.E69

The DGN method (69) always takes the descendent direction that satisfies a linear search method. It has slow convergence when the parameters to be fitted are far from the solution, but the method maximizes its advantages when it is close to the solution . For that reason, DGN should be used when the parameters are close to their optimal values.

Clearly, α defines a fraction of the step taken from GN to update the x parameters. In this work, we select α by a backtracking strategy, whereby α is reduced from unity until an acceptable value for γ(x) is found.

## 4. MATLAB algorithms and applications

In this section, algorithms for the rational approximation of frequency domain responses in MATLAB environment are provided. Rational fitting of synthetics functions through Bode, DGN, LS or Levy, IRLS, LM and VF are presented. Finally, the techniques are applied to the rational approximation of the frequency-dependent parameters corresponding to a single-phase transmission line.

### 4.1. Bode and DGN algorithms

Bode and DGN routines are presented in this section. The concept is simple, first Bode is implemented for the curve fitting process of a synthetic function, then DGN is used to improve the accuracy of the rational approximation given by Bode. A synthetic function is selected in the main program “Fitting_Bode_DGN.m” presented in Table 2. After application of the Bode routine “Bode_process.m,” shown in Table 3, a set of poles, residues and a constant term are obtained. These results are the initial values for the DGN routine “Damped_Gauss_Newton.m” presented in Table 4.

 %================================================%        Main program of BODE-DGN%   Authors: Eduardo Salvador Bañuelos Cabral%       José Alberto Gutiérrez Robles%          Bjørn Gustavsen%===============================================clcclear allclose all%% Initial SettingsNs = 400;         % Number of samplesf = logspace(-2,10,Ns); % Frequency (Hz)w = 2.*pi.*f;      % Frequency (rad/seg)s = 1j*w;         % Complex Frequency%% Synthetic functionschoice = menu('CHOOSE A FUNCTION','Function F1(s)',…   'Function F2(s)','Function F3(s)','Function F4(s)',…   'Function F5(s)','Function F6(s)');if choice == 1  Fs = (110.*s./((s+10).*(s+100)));  tol = 2.2; % Tolerance in decibels to set a new pole and/or new zero  ite = 40; % Number of iteration in DGN methodendif choice == 2   Fs = ((s+10).*(s+100))./((s+14580).*(s+10355550));   tol = 2.2; % Tolerance in decibels to set a new pole and/or new zero   ite = 40; % Number of iteration in DGN methodendif choice == 3   Fs = ((s+10655).*(s+10))./((s+148450).*(s+198545852).*(s+155222220));   tol = 2.2; % Tolerance in decibels to set a new pole and/or new zero   ite = 40; % Number of iteration in DGN methodendif choice == 4   Fs = (10124).*s.*(s+35.7).*(s+88.9)./((s+100.5).*(s+220.7).*(s+5900).*(s+1370).*(s+21000));   tol = 2; % Tolerance in decibels to set a new pole and/or new zero   ite = 50; % Number of iteration in DGN methodendif choice == 5   Fs = ((s+79).*(s+1045))./((s+1458).*(s+103555).*(s+127710355).*(s+1244103555));   tol = 1.8; % Tolerance in decibels to set a new pole and/or new zero   ite = 40; % Number of iteration in DGN methodendif choice == 6   Fs = ((s+64518).*(s+8451629).*(s+312).*(s+54841216192))./((s+456).*(s+7852).*    (s+982365).*…(s+93256888).*(s+79325684588536));   tol = 2; % Tolerance in decibels to set a new pole and/or new zero   ite = 40; % Number of iteration in DGN methodend%% Bode method[P,Z,k] = Bode_process(Fs,f,Ns,tol); % Bode subroutineas = k.*poly(Z); bs = poly(P); % Polynomials[r,p,ks] = residue(as,bs); % Poles, residues and constant termTF=isempty(ks);if(TF==1);ks=0;end%% Damped Gauss-NewtonXmcp = [r; p; ks]; % Poles, residues and constant term from BodeNp = length(p);  % Number of poles[KGN,RGN,PGN,Ff] = Damped_Gauss_Newton(Np,Ns,Fs,s,Xmcp,ite);Fs_fitDGN = zeros(1,Ns);for k = 1:length(PGN)  Fs_fitDGN = Fs_fitDGN + (RGN(k)./(s.' - PGN(k))).';endFs_fitDGN = Fs_fitDGN + KGN;error = abs(Fs - Fs_fitDGN);%% Plotsfigure(3)loglog(f,abs(Fs),'k',f,abs(Fs_fitDGN),'r‐‐',f,error,'‐‐b','LineWidth',2);xlabel('Frequency [Hz]'); ylabel('Magnitude [p.u.]');legend('Data','DGN','Deviation',2)figure (4)semilogy(Ff,'‐‐b','LineWidth',2)xlabel('Iteration count'); ylabel('RMS-error');

### Table 2.

“Fitting_Bode_DGN.m”.

 %================================================%        BODE PROCESS PROGRAM%   Authors: Eduardo Salvador Bañuelos Cabral%       José Alberto Gutiérrez Robles%          Bjørn Gustavsen%================================================% Inputs      % ‐‐‐ Fs, function to be fitted      % ‐‐‐ f, frequency (Hz)      % ‐‐‐ Ns, number of samples      % ‐‐‐ tol, tolerance% Outputs      % ‐‐‐ P, Poles      % ‐‐‐ Z, Zeros      % ‐‐‐ k, Constant termfunction [P,Z,k] = Bode_process(Fs,f,Ns,tol)P = []; Z = [];                            % Initialize P and Z like empty matricesFsdb = 20*log10(abs(Fs)); % Function in decibelsk0db = Fsdb(1);                        % k0 in decibelsk0 = 10.^(k0db./20);             % k0 in magnitudeFsdb1 = ones(1,Ns).*k0db; % Constant term% Plot of the function in decibelsfigure(1)semilogx(f,Fsdb,'k',f,Fsdb1,'r‐‐','LineWidth',2)ylabel('Decibels'), xlabel('Frequency [Hz]'), hold onlegend('Data','Bode',2)c = 1; h = 1; r = 0; % CountersFsdb2 = zeros(1,Ns); % Initialize Fsdb2while (r < Ns) % Bode process algorithm            Fsfitdb = Fsdb1 + Fsdb2;            semilogx(f,Fsfitdb,'r‐‐','LineWidth',2)            pause(0.5)            for r = 1:1:Ns                     error = abs((Fsdb(r))-(Fsfitdb(r))); % Deviation                     if error >= (tol)                                              % Tolerance                           % New zero                           if (Fsfitdb(r)< Fsdb(r)); Z(h)=f(r); h=h+1; break; end                           % New pole                           if (Fsfitdb(r)> Fsdb(r)); P(c)=f(r); c=c+1; break; end                     end           end           % Contribution of each zero in decibels           Arg1 = 0;           for kp = 1:length(Z)                    Arg1 = Arg1 + 20.*log10(abs(1+1i*f/Z(kp)));           end           % Contribution of each pole in decibels           Arg2 = 0;           for kp = 1:length(P)                   Arg2 = Arg2 - 20.*log10(abs(1+1i*f/P(kp)));           end           Fsdb2 = Arg1 + Arg2;end% Poles and zerosP = -P*2*pi;Z = -Z*2*pi;% Construct the constant termNum=1;for k=1:length(P);Num=Num*P(k);endDen=1;for k=1:length(Z);Den=Den*Z(k);endk = abs(Num)*k0/abs(Den);% Plotsfigure(2)semilogx(f,Fsdb,'k',f,Fsfitdb,'‐‐r','LineWidth',2);xlabel('Frequency [Hz]'); ylabel('Decibels')legend('Data','Bode',2)

### Table 3.

“Bode_process.m”.

 %================================================%                         DAMPED GAUSS NEWTON METHOD%           Authors: Eduardo Salvador Bañuelos Cabral%                        José Alberto Gutiérrez Robles%                                     Bjørn Gustavsen%===============================================% Inputs      % ‐‐‐ Np, number of poles      % ‐‐‐ Ns, Number of samples      % ‐‐‐ Fs, function to be fitted      % ‐‐‐ s, complex frequency (rad/s)      % ‐‐‐ X, initial poles, zeros and constant term      % ‐‐‐ ite,iterations% Outputs      % ‐‐‐ Ks, constant term      % ‐‐‐ Rs, residues      % ‐‐‐ Ps, poles      % ‐‐‐ Ff, objective function (deviation)function [Ks,Rs,Ps,Ff] = Damped_Gauss_Newton(Np,Ns,Fs,s,X,ite)a = ones(Ns,1);                              % Vector columnJ1 = zeros(Ns,Np);                       % Matrix size (J1)J2 = zeros(Ns,Np);                       % Matrix size (J2)Jn = zeros(2*Ns,2*Np+1);         % Matrix size (Jn)At = zeros(2*Np+1,2*Np+1);    % Matrix size (At)en = zeros(2*Ns,1);                    % Vector size (en)epn = zeros(2*Ns,1);                  % Vector size (epn)Ff = zeros(ite,1);                       % Vector size (Ff)Euclidian = zeros(1,1+2*Np); % Vector sizefor ki = 1:ite % Damped Gauss Newton methodology        R = X(1:Np); % Residues        P = X(Np+1:2*Np); % Poles        K = X(2*Np+1); % Constant term        Fa = zeros(1,Ns); % Set the approximation to the function        for k = 1:length(P)                Fa = Fa + (R(k)./(s.' - P(k))).';        end        Fa = Fa + K; % Construct the approximation to the function        error = (Fa.' - Fs.'); % Deviation        for n = 1:1:Np % Loop to construct the Jacobian                J1(1:Ns,n) = 1./(s-P(n));                J2(1:Ns,n) = R(n)./((s-P(n)).^2);        end        J = [J1 J2 a];                     % Jacobian        [Xmax Ymax] = size(J); % Matrix size (J)        Jr = real(J);                     % Real part of vector J        Ji = imag(J);                     % Imaginary part of vector J        er = real(error);            % Real part of vector error        ei = imag(error);            % Imaginary part of vector error        km = 1; % Counter        for k = 2:2:2*Xmax % Interleaved                Jn(k-1,:) = Jr(km,:);                Jn(k,:) = Ji(km,:);                en(k-1,1) = er(km);                en(k,1) = ei(km);                km = km+1;        end        F = ((norm(en,2)^2)); % Objective function        Ff(ki,1) = F;                    % Storage the objective function must tend to zero       [Q,R] = qr(Jn);                % Matrix Q and R of Jn        Jn = R;                                 % It updates matrix Jn        Gf = (Jn.'*Q.')*en;      % Gradient (QR decomposition)        Hess = Jn.'*Jn;               % Hessian (QR decomposition)        for col = 1:Ymax % Euclidian norm                Euclidian(col) = norm(Hess(:,col),2);                At(:,col) = Hess(:,col)./Euclidian(col);        end        h = (At)\-Gf;               % Solution for the system (Ax = b)        h = h./Euclidian.'; % Real solution        stop = 1;                        % Variable to stop the next loop        al = 1;                             % Variable to weigh the approximation        while (stop == 1)                Xp = X + al*h;                % New coefficients (without updating x)                Rp = Xp(1:Np);              % Residues                Pp = Xp(Np+1:2*Np);  % Poles                Kp = Xp(2*Np+1);         % Constant term                Fap = zeros(1,Ns);    % Set the new approximation of the function                for k = 1:length(P)                        Fap = Fap + (Rp(k)./(s.' - Pp(k))).';                end                Fap = Fap + Kp;            %New approximation of the function                ep = (Fap.' - Fs.'); % Deviation                epr = real(ep);          % Real part of vector ep                epi = imag(ep);          % Imaginary part of vector ep                km = 1; % Counter                for k = 2:2:2*Xmax % Interleaved                         epn(k-1,1) = epr(km);                         epn(k,1) = epi(km);                         km = km+1;                 end                 Fp = ((norm(epn,2)^2)); % Objective function                 % Updating process                 Erel = 1e-4*al*h.'*Gf; % Relative error to stop the while process                 if (Fp < F + Erel)                        X = X + al*h; % Final approximation                        stop = 0; % Go out the while                 else                         al = al*0.9; % Weigh to update X                         if (al < 1e-30)                                stop = 0;                          else                                stop = 1;                          end                 end           endendRs = X(1:Np);            % ResiduesPs = X(Np+1:2*Np); % PolesKs = X(2*Np+1);        % Constant term

### Table 4.

“Damped_Gauss_Newton.m”.

In Figure 5 the asymptotic behavior of the Bode routine is shown, and in Figure 6 the final rational approximation given by Bode is presented. Following, DGN routine improves the curve fitting process, this is presented in Figure 7. Moreover, the RMS-error on each iteration is shown in Figure 8.

### 4.2. LS or Levy algorithm

In this section the LS method or Levy method routine is introduced. In the main program “Fitting_OLS.m,” Table 5, a synthetic function can be selected by the user. Then the LS routine “Least_Squares_Method.m,” Table 6, is implemented for the rational approximation of this function.

 %====================================================%                 Main program of ORDINARY LEAST SQUARES%              Authors: Eduardo Salvador Bañuelos Cabral%                           José Alberto Gutiérrez Robles%                                     Bjørn Gustavsen%====================================================clcclear allclose all%% Initial SettingsNs = 400;                               % Number of samplesf = logspace(-2,6,Ns); % Frequency "Hertz"w = 2.*pi.*f;                      % Frequency "rad/seg"s = 1i*w;                               % Complex Frequency%% Synthetic functionschoice = menu('CHOOSE A FUNCTION','Function F1(s)','Function F2(s)','Function F3(s)');if choice == 1Np = 2;Fs = 0.01+(210*(s)./((s+10).*(s+100)));endif choice == 2      Np = 5;      Fs = (-81+18e6*s+27e6*s.^2+62e3*s.^3+42*s.^4+5e-7*s.^5)./…               (1+54e5*s+36e7*s.^2+83e4*s.^3+90*s.^4+36e-3*s.^5);endif choice == 3      Np = 6;      Fs =(100./(s+200))+(250./(s+2000))+(1500./(s+(10-1j*100)))+(1500./(s+(10+1j*100)))+…              (800./(s+(1000+1j*30000)))+(800./(s+(1000-1j*30000)));end%% Ordinary Least Squares[Ks,Rs,Ps,x,A1,an,bn] = Least_Squares_Method(Np,Ns,Fs,s);% Fitting using partial fractionsFs_fit = zeros(1,Ns);for i = 1:Np        Fs_fit = Fs_fit + (Rs(i)./(s.' - Ps(i))).';endFs_fit = Fs_fit + Ks;error = abs((Fs.' - Fs_fit.' ));%% Plotsfigure(1)loglog(f,abs(Fs),'-k',f,abs(Fs_fit),'‐‐r',f,error,'-.b','LineWidth',2);legend('Data','Ordinary Least Squares','Deviation',3);xlabel('Frequency [Hz]'); ylabel('Magnitude [p.u.]')

### Table 5.

“Fitting_OLS.m”.

 %==========================================================%                  Function to implement ORDINARY LEAST SQUARES%                     Authors: Eduardo Salvador Bañuelos Cabral%                               José Alberto Gutiérrez Robles%                                              Bjørn Gustavsen%==========================================================% Inputs      % ‐‐‐ Np, number of poles      % ‐‐‐ Ns, number of samples      % ‐‐‐ Fs, function to be fitted      % ‐‐‐ s, complex frequency (rad/s)% Outputs      % ‐‐‐ Ks, constant term      % ‐‐‐ Rs, residues      % ‐‐‐ Ps, poles      % ‐‐‐ x, vector of coefficients (polynomials)      % ‐‐‐ A1, matrix to evaluate the fitting      % ‐‐‐ an, coefficients of the numerator      % ‐‐‐ bn, coefficients of the denominatorfunction [Ks,Rs,Ps,x,A1,an,bn] = Least_Squares_Method(Np,Ns,Fs,s)a = ones(Ns,1);                                % Column vectora1 = zeros(Ns,Np);                         % Matrix size (s^n)a2 = zeros(Ns,Np);                         % Matrix size (Fs*s^n)An = zeros(2*Ns,1+2*Np);            % Matrix size (real and imaginary part)Bn = zeros(2*Ns,1);                        % Vector size (real and imaginary part)Euclidian = zeros(1,1+2*Np);  % Vector size% Construction of vector b.b = Fs.';% Construction of matrix A.for n = 1:1:Np        a1(1:Ns,n) = s.^((Np+1)-n);        a2(1:Ns,n) = -b.*a1(1:Ns,n);endA1 = [a1 a];    % Matrix to evaluate the fittingA = [a1 a a2]; % Matriz of the system[Xmax Ymax] = size(A); % Size for matrix AAr = real(A);                    % Real part of matrix AAi = imag(A);                    % Imaginary part of matrix Abr = real(b);                    % Real part of vector Bbi = imag(b);                    % Imaginary part of vector Bkm = 1; % Construction of matrix An and vector Bn (interleaved)for k = 2:2:2*Xmax        An(k-1,:) = Ar(km,:);        An(k,:) = Ai(km,:);        Bn(k-1,1) = br(km);        Bn(k,1) = bi(km);        km = km+1;end[Q,R] = qr(An,0); % QR decomposition - Matrix Q and R of AnAt = R;                        % It updates the matrix AnB = Q.'*Bn;               % It updates the vector Bnfor col = 1:Ymax     Euclidian(col) = norm(At(:,col),2); % Euclidian norm     At(:,col) = At(:,col)./Euclidian(col);endx = At\B;                        % Solution for the system (Ax = b)x = x./Euclidian.'; % Real solutionan = x(1:Np+1);                     % Numerator coefficientsbn = [x(Np+2:2*Np+1);1]; % Denominator coefficients% Poles, residues and constant term[Rs,Ps,Ks] = residue(an,bn);TF = isempty(Ks);if (TF==1)       Ks = 0;end

### Table 6.

“Least_Squares_Method.m”.

Figure 9 shows the synthetic frequency response data and the approximation given by the LS routine, additionally, the deviation in terms of absolute error is presented. Figure 9.Synthetic frequency response data and the approximation given by the LS routine.

### 4.3. IRLS algorithm

Because in practice the Noda method and the SK method correspond to a specific weighting function in the IRLS technique, only the latter is presented. In the beginning of the implementation a synthetic function is selected in the main program “Fitting_IRLS.m,” presented in Table 7. Afterwards, the algorithm of the method “IRLS.m,” Table 8, is implemented for the rational approximation of the selected function.

 %=======================================================%       Main program of ITERATIVELY REWEIGHTED LEAST SQUARES%                   Authors: Eduardo Salvador Bañuelos Cabral%                                José Alberto Gutiérrez Robles%                                           Bjørn Gustavsen%=======================================================clcclear allclose all%% Initial SettingsNs = 500;                               % Number of samplesf = logspace(-2,6,Ns); % Frequency (Hz)w = 2.*pi.*f;                      % Frequency (rad/seg)s = 1i*w;                               % Complex Frequency%% Synthetic functionschoice = menu('CHOOSE A FUNCTION','Function F1(s)',…         'Function F2(s)','Function F3(s)','Function F4(s)',…         'Function F5(s)','Function F6(s)');if choice == 1      Np = 12;      Fs =(0.83*s.^12-0.35*s.^11+0.32*s.^10+60000*s.^9+600.36*s.^8+650.23*s.^7+…               54.5*s.^6+45.6*s.^5 +0.1*s.^4-1240*s.^3+80025*s.^2+5547.*s+5008)./…              (s.^12-0.89*s.^11-0.23*s.^10-4565.*s.^9+34569.02*s.^8-55450*s.^7+…              60*s.^6+1005*s.^5+8050*s.^4+7502*s.^3+10*s.^2+1014.*s-50);endif choice == 2      Np = 9;      Fs = (6000*s.^9+600*s.^8+650*s.^7+54*s.^6+457*s.^5+100*s.^4+1240*s.^3+8000*s.^2+…               55400.*s+5000)./(450*s.^9+3500*s.^8+500*s.^7-600*s.^6+1000*s.^5+700*s.^4-…               7500*s.^3+1000*s.^2+500.*s);endif choice == 3      Np = 7;Fs =(514*s.^7-364.25*s.^6-0.35*s.^5+635*s.^4-20802*s.^3+5304*s.^2+4520.*s-…        18020.025 )./(241*s.^7-32.02*s.^6+538.23*s.^5-2588*s.^4-22560*s.^3-…        604*s.^2+4150.*(s)-21052.31 );endif choice == 4      Np = 3;      Ps = [-0.52; -0.12; -0.04];      Rs = [-2.18; -7192.25; 20.58];      Ks = 5;      Fs = zeros(1,Ns);for k = 1:length(Ps)        Fs = Fs + (Rs(k)./(s.' - Ps(k))).';endFs = Fs + Ks;endif choice == 5      Np = 11;      Ps = 1.0e+004.*[-.88; -5.42; - 27.28; 0.88; -9665; 235.7; …                                         -0.26; -3.87; 46.32; -0.13; -3756 ];      Rs = 1.0e+005.*[834; 22593; 893; 2653; 654; 32; …                                        44.14; 6405; 136.79; 0.12; 125];      Ks = 0;      Fs = zeros(1,Ns);      for k = 1:length(Ps)              Fs = Fs + (Rs(k)./(s.' - Ps(k))).';      end      Fs = Fs + Ks;endif choice == 6      Np = 18;      Ps = [-4500; -41000; -100+1i*5000; -100-1i*5000; -120+1i*15000; -120-1i*15000; …              -3000+1i*35000; -3000-1i*35000; -200+1i*45000; -200-1i*45000; -1500+1i*45000;…              -1500-1i*45000; -500+1i*70000; -500-1i*70000; -1000+1i*73000; -1000-1i*73000;…              -2000+1i*90000; -2000-1i*90000];      Rs = [-3000; -83000; -5+1i*7000; -5-1i*7000; -20+1i*18000; -20-1i*18000;…              6000 + 1i*45000; 6000-1i*45000; 40+1i*60000; 40-1i*60000; 90+1i*10000;…              90-1i*10000; 50000+1i*80000; 50000-1i*80000; 1000+1i*45000; 1000-1i*45000;…             -5000+1i*92000; -5000-1i*92000];     Ks = 0.2;     Fs = zeros(1,Ns);     for k = 1:length(Ps)             Fs = Fs + (Rs(k)./(s.' - Ps(k))).';endFs = Fs + Ks;end%% Iteratively Reweighted Least Squaresfw = 1; % weighting function (1-6)[Ks,Rs,Ps,x,A1,an,bn,k1,Err] = IRLS(Np,Ns,Fs,s,fw);% Fitting using partial fractionsFs_fit = zeros(1,Ns);for i = 1:Np        Fs_fit = Fs_fit + (Rs(i)./(s.' - Ps(i))).';endFs_fit = Fs_fit + Ks ;e = abs((Fs.' - Fs_fit.')); % Deviation of the fitting process%% Plotsfigure(1)semilogy(Err(1:length(Err)),'-b','LineWidth',2)xlabel('Itetation count');ylabel('RMS-error')figure(2)loglog(f,abs(Fs),'-k',f,abs(Fs_fit),'‐‐r',f,e,'‐‐b','LineWidth',2),legend('Data','IRLS','Deviation')xlabel('Frequency [Hz]'); ylabel('Magnitude [p.u.]');

### Table 7.

“Fitting_IRLS.m”.

 %=====================================================%               ITERATIVELY REWEIGHTED LEAST SQUARES METHOD%                 Authors: Eduardo Salvador Bañuelos Cabral%                          José Alberto Gutiérrez Robles%                                      Bjørn Gustavsen%=====================================================% Inputs      % ‐‐‐ Np, number of poles      % ‐‐‐ Ns, number of samples      % ‐‐‐ Fs, function to be fitted      % ‐‐‐ s, complex frequency (rad/s)      % ‐‐‐ fw, weighting function% Outputs      % ‐‐‐ Ks, constant term      % ‐‐‐ Rs, residues      % ‐‐‐ Ps, poles      % ‐‐‐ x, coefficients vector (polynomials)      % ‐‐‐ A1, matrix to evaluate the fitting      % ‐‐‐ an, numerator coefficients      % ‐‐‐ bn, denominator coefficients      % ‐‐‐ k1, minimum error      % ‐‐‐ Error, RMS-errorfunction [Ks,Rs,Ps,x,A1,an,bn,k1,Error] = IRLS(Np,Ns,Fs,s,fw)a = ones(Ns,1);                               % Vector columna1 = zeros(Ns,Np);                        % Matrix size (s^n)a2 = zeros(Ns,Np);                        % Matrix size (Fs*s^n)An = zeros(2*Ns,1+2*Np);          % Matrix size considering real and imaginary partBn = zeros(2*Ns,1);                     % Vector size considering real and imaginary partEuclidian = zeros(1,1+2*Np); % Vector sizeW = eye(length(Ns*2));               % Weighting matrixb = Fs.'; % Construction of vector b.for n = 1:1:Np % Construction of matrix A.        a1(1:Ns,n) = s.^((Np+1)-n);        a2(1:Ns,n) = -b.*a1(1:Ns,n);endA1 = [a1 a];     % Matrix to evaluate the fittingA = [a1 a a2]; % Matriz of the system[Xmax Ymax] = size(A); % Size for matrix AAr = real(A);                    % Real part of matrix AAi = imag(A);                    % Imaginary part of matrix Abr = real(b);                    % Real part of vector Bbi = imag(b);                    % Imaginary part of vector B% Construction of matrix An and vector Bn (interleaved)km = 1;for k = 2:2:2*Xmax        An(k-1,:) = Ar(km,:);        An(k,:) = Ai(km,:);        Bn(k-1,1) = br(km);        Bn(k,1) = bi(km);        km = km+1;end[Q,R] = qr(An,0);% Matrix Q and R of AnAt = R;                       % It updates matrix AnB = Q.'*Bn;              % It updates vector Bnfor col = 1:Ymax    % Applying the Euclidian norm to At     Euclidian(col) = norm(At(:,col),2);      At(:,col) = At(:,col)./Euclidian(col);endxn(:,1) = At\B;                                     % Solution for the system (Ax = B)xn(:,1) = xn(:,1)./Euclidian.'; % Real solution% weighting iterativeError(1) = 1e2; % Initial errorDisc = 1e2; % Discriminatingn=1; % Counter to storage de errork0=1; % Counter of iterationswhile Disc < 0 || Error(n) > 1e-3 || k0 < 13          an = xn(1:Np+1,k0);                      % Numerator coefficients          bn = [xn(Np+2:2*Np+1,k0);1]; % Denominator coefficients          Fsfit = (A1*an)./(A1*bn);        % Fitting evaluation          vres = (Fsfit - b);                        % Direct error          n=n+1;                                                 % Counter to storage de error          Error(n)=sum(abs(vres));        % Error          Disc = Error(n)-Error(n-1);   % Discriminating          vres_r = real(vres);                  % Real part of the vector vres          vres_i = imag(vres);                  % Imaginary part of the vector vres          km = 1; %(interleaved          for k = 2:2:2*Xmax                  res(k-1,1) = vres_r(km);                  res(k,1) = vres_i(km);                  km = km+1;          end          Ds = std(res); % Standard deviation          res = res./Ds; % Vector "res" with standard deviation          % WEIGHING FUNCTIONS          if fw ==1, w = diag(W).*abs(res) + 1e-130; end          if fw ==2, w = diag(W).*(sqrt(1+((res.^2)))-1) + 1e-130; end          if fw ==3, w = diag(W).*1./(abs(res).^(1.2-2)) + 1e-130; end          if fw ==4, w = diag(W).*((res).^2)./(1+(res).^2) + 1e-130; end          if fw ==5, w = diag(W).*((res).^2)./sqrt(1+(res).^2)+ 1e-130; end          if fw ==6, w = diag(W).*(abs(res).^(1.2))./1.2 + 1e-130; end          W = diag(w);                 % It updates the matrix W whit weighting functions          Bn_p = W*Bn;                 % Weighting An          An_p = W*An;                 % Weighting Bn          [Q,R] = qr(An_p,0); % Matrix Q and R of An_p          At = R;                             % It updates the matrix An_p          B = Q.'*Bn_p;               % It updates the matrix Bn_p          for col = 1:Ymax % Euclidian norm                 Euclidian(col) = norm(At(:,col),2);                 At(:,col) = At(:,col)./Euclidian(col);          end          k0=k0+1;                                                       % Number of iterations          xn(:,k0) = At\B;                                       % Solution for the system (Ax = b)          xn(:,k0) = xn(:,k0)./Euclidian.'; % Real solution          if k0 > 20, break, end % End the while, if k0 > 20end[k1,k2]=min(Error); % Position of the minimum errorx=xn(:,k2);                    % Coefficients with minimum error% Poles, residues and constant terman = x(1:Np+1);                     % Numerator coefficientsbn = [x(Np+2:2*Np+1);1]; % Denominator coefficients[Rs,Ps,Ks] = residue(an,bn);TF = isempty(Ks);if (TF==1)       Ks = 0;end

### Table 8.

“IRLS.m”.

In Figure 10 the synthetic frequency response behavior and the rational approximation given by IRLS is presented, together with its deviation in terms of absolute error. Moreover, the RMS-error on each iteration is shown in Figure 11. Figure 10.Rational approximation given by the IRLS method for the synthetic function. Figure 11.RMS-error on each iteration given by the IRLS method.

### 4.4. LM algorithm

In this section the LM routine is introduced. First, a synthetic function can be selected in the main program “Fitting_LM_Polynomials.m,” presented in Table 9. Afterwards, the algorithm of the method “LM_method.m,” Table 10, is implemented for the rational approximation of the selected function.

 %=======================================================%        Main program of LEVENBERG-MARQUARDT (Polynomials)%                  Authors: Eduardo Salvador Bañuelos Cabral%                             José Alberto Gutiérrez Robles%                                           Bjørn Gustavsen%=======================================================clcclear allclose all%% Initial SettingsNs = 500;                               % Number of samplesf = logspace(-2,8,Ns); % Frequency "Hertz"w = 2.*pi.*f;                      % Frequency "rad/seg"s = 1i*w;                               % Complex Frequency%% Synthetic functionschoice = menu('CHOOSE A FUNCTION','Function F1(s)','Function F2(s)',…'Function F3(s)','Function F4(s)','Function F5(s)');if choice == 1      ite = 150; % Number of iterations      Np = 7; % Order of the proposed function      Fs = (-364*s.^6-s.^5+635*s.^4-20802*s.^3+5304*s.^2+4520.*s-18020)./…               (241*s.^7-32*s.^6+538*s.^5-2588*s.^4-22560*s.^3-604*s.^2+4150.*s-21052);endif choice == 2      ite = 50; % Number of iterations      Np = 12; % Order of the proposed function      Fs = (s.^12-0.4*s.^11+0.3*s.^10+60003*s.^9+600*s.^8+650*s.^7+54*s.^6+45*s.^5+…               0.1*s.^4-1240*s.^3+80025*s.^2+5547*s+5008)./(s.^12-0.9*s.^11-0.2*s.^10-…               4565*s.^9+34569*s.^8-55450*s.^7+60*s.^6+1005*s.^5+8050*s.^4+7502*s.^3+…               10*s.^2+1014*s-50);endif choice == 3      ite = 50; % Number of iterations      Np = 14; % Order of the proposed function      Fs = 2.*pi*(652*s.^14-0.8*s.^13-0.8*s.^12-0.4*s.^11+25487*s.^10+600003*s.^9+…              6040*s.^8-60050*s.^7+54*s.^6+45*s.^5+0.1*s.^4-1240*s.^3+8025*s.^2+5547*s+508)./…              (2365*s.^14-8362*s.^13-8*s.^12-0.9*s.^11-65547*s.^10-4565*s.^9+34569*s.^8-…              55450*s.^7+60*s.^6-1005*s.^5+850*s.^4-7502*s.^3-10*s.^2-1014*s+75230);endif choice == 4      ite = 50; % Number of iterations      Np = 3; % Order of the proposed function      Ps = [-0.5; -0.1; -0.05];      Rs = [-2.2; -7192; 20];      % Process to obtain Fs from poles and residues      Ks = 5;      Fs = zeros(1,Ns);      for k = 1:length(Ps)              Fs = Fs + (Rs(k)./(s.' - Ps(k))).';      end      Fs = Fs + Ks;endif choice == 5      ite = 100; % Number of iterations      Np = 11; % Order of the proposed function      Ps = 1.0e4*[-0.9;-5.4;-27.3;0.9;-9665;235;-0.3;-3.9;46.3;-0.1;-3756];      Rs = 1.0e5*[834;22593;893;2653;654;32;44;6405;136;0.1;125];      % Process to obtain Fs from poles and residues      Ks = 0;      Fs = zeros(1,Ns);      for k = 1:length(Ps)              Fs = Fs + (Rs(k)./(s.' - Ps(k))).';      end      Fs = Fs + Ks;end%% Levenberg Marquardt[Ks,Rs,Ps,x,A1,an,bn,Ff] = LM_method(Np,Ns,Fs,s,ite);% Fitting using partial fractionsFs_fit = zeros(1,Ns);for i = 1:Np        Fs_fit = Fs_fit + (Rs(i)./(s.' - Ps(i))).';endFs_fit = Fs_fit + Ks ;% Deviation in the fitted processerror = abs((Fs.' - Fs_fit.' ));%% Plotsfigure(1)loglog(f,abs(Fs),'-k',f,abs(Fs_fit),'-.r',f,error,'-.b','LineWidth',2);legend('Data','LM (Polynomials)','Deviation');xlabel('Frequency [Hz]'); ylabel('Magnitude [p.u.]');figure(2)semilogy(Ff,'‐‐b','LineWidth',2);xlabel('Iteration count');ylabel('RMS-error');

### Table 9.

“Fitting_LM_Polynomials.m”.

 %==========================================================%      Function to implement LEVENBERG-MARQUARDT (Polynomials)%                    Authors: Eduardo Salvador Bañuelos Cabral%                                José Alberto Gutiérrez Robles%                                            Bjørn Gustavsen%==========================================================% Inputs      % ‐‐‐ Np, number of poles      % ‐‐‐ Ns, number of samples      % ‐‐‐ Fs, function to be fitted      % ‐‐‐ s, complex frequency (rad/s)      % ‐‐‐ ite, number of iterations% Outputs      % ‐‐‐ Ks, constant term      % ‐‐‐ Rs, residues      % ‐‐‐ Ps, poles      % ‐‐‐ x, vector of coefficients (polynomials)      % ‐‐‐ A1, matrix to evaluate the fitting      % ‐‐‐ an, coefficients of the numerator      % ‐‐‐ bn, coefficients of the denominator      % ‐‐‐ Ff, objective function (deviation)function [Ks,Rs,Ps,x,A1,an,bn,Ff] = LM_method(Np,Ns,Fs,s,ite)x = ones((2*Np)+1,1);                 % Initial valuesa = ones(Ns,1);                              % Vector columna1 = zeros(Ns,Np);                       % Matrix size (s^n)J1 = zeros(Ns,Np+1);                   % Matrix size (J1)J2 = zeros(Ns,Np);                       % Matrix size (J1)At = zeros(2*Np+1,2*Np+1);     % Matrix size (At)en = zeros(2*Ns,1);                     % Vector size (en)epn = zeros(2*Ns,1);                   % Vector size (enp)Ff = zeros(ite,1);                        % Vector size (Ff)Euclidian = zeros(1,1+2*Np); % Euclidian norm% Construction of the matrix A1for n = 1:1:Np        a1(1:Ns,n) = s.^((Np+1)-n);endA1 = [a1 a];v = 2; % LM coefficient (m = m*v)u = 5; % LM coefficient (m = m/u)% Levenberg-Marquardt methodfor ki = 1:ite % Main loop        x1 = x(1:Np+1); % Coefficients of the numerator        x2 = [x(Np+2:2*Np+1);1]; % Coefficients of the denominator        Yg = (A1*x1)./(A1*x2); % Fitting evaluation        error = (Yg - Fs.'); % Deviation        % Process to obtain the Jacobian (J)        for k = 1:1:Np               J1(:,k) = (s.^((Np+1)-k))./((A1*x2).');               J2(:,k) = (-(A1*x1).'.*s.^((Np+1)-k))./((A1*x2).^2).';        end        J1(:,Np+1) = 1./((A1*x2).');        J = [J1 J2];        [Xmax Ymax] = size(J); % Matrix size (J)        Jr = real(J);                    % Real part of vector J        Ji = imag(J);                    % Imaginary part of vector J        er = real(error);           % Real part of vector e        ei = imag(error);           % Imaginary part of vector e        km = 1; % Interleaved        for k = 2:2:2*Xmax               Jn(k-1,:) = Jr(km,:);              Jn(k,:) = Ji(km,:);              en(k-1,1) = er(km);              en(k,1) = ei(km);              km = km+1;        end        F = ((norm(en,2)^2)); % Objective function tends to zero to converge        Ff(ki,1) = F; % Storage the RMS-error (objective function)        [Q,R] = qr(Jn);          % Matrix Q and R of Jn        Jn = R;                            % It updates matrix Jn        Gf = (Jn.'*Q.')*en; % Gradient (QR decomposition)        Hess = Jn.'*Jn;           % Hessian (QR decomposition)        if ki == 1                        % Inicial value for m        m = max(diag(Hess));        end        paro = 0; % Variable to enter into the while        while (paro==0)                 I = diag(diag(Hess)); % Construction of matrix I, diagonal of Hessian                 Hess_mod = (Hess + m*I); % Modified Hessian                 for col = 1:Ymax % Loop to normalize At                         Euclidian(col) = norm(Hess_mod(:,col),2); % Euclidian norm                          At(:,col) = Hess_mod(:,col)./Euclidian(col);                 end                 h = (At)\-Gf;                % Solution for the system (Ax = b)                 h = h./Euclidian.';  % Real solution                 xnew = x + h;                  % New coefficients (without updating x)                 % == Updated coefficients are evaluated, if it meets the assessment then x is updated == %                 x1new = xnew(1:Np+1);                         % Coefficients of the numerator                 x2new = [xnew(Np+2:2*Np+1);1];    % Coefficients of the denominator                 Ygnew = (A1*x1new)./(A1*x2new);  % Fitting evaluation                 ep = (Ygnew - Fs.');                              % Deviation                 epr = real(ep);                                       % Real part of vector ep                 epi = imag(ep);                                       % Imaginary part of vector ep                 km = 1; % Interleaved                 for k = 2:2:2*Xmax                         epn(k-1,1) = epr(km);                         epn(k,1) = epi(km);                         km = km+1;                 end                 Fp = ((norm(epn,2)^2)); % Objective function                 % Updating of m                 if (F < Fp)                        m = m*v;                 end                                  if (F >= Fp) % Condition to go out the while                        x = x + h;                        m = m/u;                        break;                 end         endendan = x(1:Np+1); % Coefficients of the numeratorbn = [x(Np+2:2*Np+1);1]; % Coefficients of the denominator% Poles, residues and constant term[Rs,Ps,Ks] = residue(an,bn);TF = isempty(Ks);if (TF==1)       Ks = 0;end

### Table 10.

“LM_method.m”.

In Figure 12 the synthetic frequency response behavior and the rational approximation given by IRLS is presented, together with its deviation in terms of absolute error. Moreover, the RMS-error on each iteration is shown in Figure 13.

### 4.5. VF algorithm

VF has become one of the main methodologies for the rational approximation of frequency domain responses, the MATLAB routine is available online , “vectfit3.m.” Here, we present an example of its application to the rational approximation of a synthetic function.

In the beginning of the implementation a synthetic function is selected in the main program “Fitting_VF.m,” presented in Table 11. Afterwards, initial poles for the fitting must be established “InitialPoles.m,” Table 12, for the VF routine. Finally in “VF.m,” VF is used (via vectfit3.m) to perform the rational approximation of the selected function. In VF.m routine, Table 13, the parameter opts.relax is set to 1. This implies that “relaxation” is being used.

 %=====================================================================%                                               Main program of VECTOR FITTING%      (vectfit3 available in:https://www.sintef.no/projectweb/vectfit/%                                                    Author: Bjørn Gustavsen)%                                Authors: Eduardo Salvador Bañuelos Cabral%                                            José Alberto Gutiérrez Robles%                                                           Bjørn Gustavsen%====================================================================clcclear allclose all%% Initial SettingsNs = 500;                               % Number of samplesf = logspace(-2,6,Ns); % Frequency (Hz)w = 2.*pi.*f;                      % Frequency (rad/seg)s = 1i*w;                               % Complex Frequency%% Synthetic functionschoice = menu('CHOOSE A FUNCTION','Function F1(s)',…          'Function F2(s)','Function F3(s)','Function F4(s)');if choice == 1      Np = 12;      Fs =(0.83*s.^12-0.35*s.^11+0.32*s.^10+60000*s.^9+600.36*s.^8+650.23*s.^7+…               54.5*s.^6+45.6*s.^5 +0.1*s.^4-1240*s.^3+80025*s.^2+5547.*s+5008)./…               (s.^12-0.89*s.^11-0.23*s.^10-4565.*s.^9+34569.02*s.^8-55450*s.^7+60*s.^6+…              1005*s.^5+8050*s.^4+7502*s.^3+10*s.^2+1014.*s-50);endif choice == 2      Np = 9;      Fs = (6000*s.^9+600*s.^8+650*s.^7+54*s.^6+457*s.^5+100*s.^4+1240*s.^3+8000*s.^2+…                55400.*s+5000 )./(450*s.^9+3500*s.^8+500*s.^7-600*s.^6+1000*s.^5+…                700*s.^4-7500*s.^3+1000*s.^2+500.*s);endif choice == 3      Np = 7;      Fs =(514*s.^7-364.25*s.^6-0.35*s.^5+635*s.^4-20802*s.^3+5304*s.^2+4520.*s-…              18020.025 )./(241*s.^7-32.02*s.^6+538.23*s.^5-2588*s.^4-22560*s.^3-…              604*s.^2+4150.*(s)-21052.31 );endif choice == 4      Np = 18;      Ps = 2*pi*[-4500; -41000; -100+1i*5000; -100-1i*5000; -120+1i*15000; …              -120-1i*15000; -3000+1i*35000; -3000-1i*35000; -200+1i*45000; -200-1i*45000;…              -1500+1i*45000; -1500-1i*45000;-500+1i*70000; -500-1i*70000; -1000+1i*73000;…              -1000-1i*73000; -2000+1i*90000; -2000-1i*90000];      Rs = [-3000; -83000; -5+1i*7000; -5-1i*7000; -20+1i*18000; -20-1i*18000;…              6000+1i*45000; 6000-1i*45000; 40+1i*60000; 40-1i*60000; 90+1i*10000; …              90-1i*10000; 50000+1i*80000; 50000-1i*80000; 1000+1i*45000; 1000-1i*45000;…              -5000+1i*92000; -5000-1i*92000];      Ks = 0.2;      Fs = zeros(1,Ns);      for k = 1:length(Ps)              Fs = Fs + (Rs(k)./(s.' - Ps(k))).';      end      Fs = Fs + Ks ;end%% Vector Fitting Method[Ps]=InitialPoles(f,Np); %Initial poles subroutine[P,C,D,E,fVF] = VF(Ps,Ns,Fs,s,20,3); % Vector Fitting MethodFsfit = zeros(1,Ns);for k = 1:length(P)        Fsfit = Fsfit + (C(k)./(s.' - P(k))).';endFsfit = Fsfit + D+ E.*s; % Fitting resulteVF = abs(Fs - Fsfit); % Deviationfigure(1)loglog(f,abs(Fs),'-k',f,abs(Fsfit),'‐‐r',f,abs(eVF),'b','LineWidth',2);xlabel('Frequency [Hz]'); ylabel('Magnitude [p.u.]')legend('Data','VF','Deviation')

### Table 11.

“Fitting_VF.m”.

 %==========================================================%        Function to implement INITIAL POLES for Vector Fitting%                   Authors: Eduardo Salvador Bañuelos Cabral%                              José Alberto Gutiérrez Robles%                                            Bjørn Gustavsen%==========================================================% Inputs       % ‐‐‐ f, frequency       % ‐‐‐ Npol, number of poles% Outputs       % ‐‐‐ Ps, polesfunction [Ps] = InitialPoles(f,Npol)% Set the initial poles (even, odd)even = fix(Npol/2);p_odd = Npol/2 - even;disc = p_odd ~= 0;% Set a real pole in case of disc == 1if disc == 0 % Even initial poles      pols = [];else                % Odd initial poles      pols = [(max(f)-min(f))/2];end% Initial complex polesbet = linspace(min(f),max(f),even);for n=1:length(bet)         alf=-bet(n)*1e-2;         pols=[pols (alf-1j*bet(n)) (alf+1j*bet(n)) ];endPs = pols.'; % Column vector of initial poles

### Table 12.

“InitialPoles.m”.

 %=====================================================================%                                           Function to use VECTOR FITTING%   (vectfit3.m available in:https://www.sintef.no/projectweb/vectfit/%                                                Author: Bjørn Gustavsen)%                             Authors: Eduardo Salvador Bañuelos Cabral%                                           José Alberto Gutiérrez Robles%                                                          Bjørn Gustavsen%=====================================================================% Inputs      % ‐‐‐ Ps, initial poles      % ‐‐‐ Ns, number of samples      % ‐‐‐ Fs, function to be fitted      % ‐‐‐ s, complex frequency (rad/s)      % ‐‐‐ ite, iterations      % ‐‐‐ Ka, kind of fitting% Outputs      % ‐‐‐ P, poles      % ‐‐‐ C, residues      % ‐‐‐ D, constant term      % ‐‐‐ E.      % ‐‐‐ RMS.function [P,C,D,E,RMS] = VF(Ps,Ns,Fs,s,Ni,Ka)weight=ones(1,Ns);  % Vector of weightsopts.relax=1;              % Use vector fitting with relaxed non-triviality constraintopts.stable=0;           % Enforce stable polesopts.asymp=Ka;           % Include both D, E in fittingopts.skip_pole=0;    % Do NOT skip pole identificationopts.skip_res=0;      % DO skip identification of residues (C,D,E)opts.cmplx_ss=1;      % Create real-only state space modelopts.spy1=0;               % No plotting for first stage of vector fittingopts.spy2=0;               % Create magnitude plot for fitting of f(s)opts.logx=1;               % Use linear abscissa axisopts.logy=1;               % Use logarithmic ordinate axisopts.errplot=1;        % Include deviation in magnitude plotopts.phaseplot=0;   % Do NOT produce plot of phase angleopts.legend=1;          % Include legends in plotsfor k = 1:Ni % Loop to make N-iterations       [SER,Ps,rmserr,fit] = vectfit3(Fs,s,Ps,weight,opts);        RMS(k) = rmserr;endP = Ps; % PolesC = SER.C; % ResiduesD = SER.D; % Constant termE = SER.E; % Proportional term

### Table 13.

“VF.m”.

Figure 14 shows the synthetic frequency response behavior and the rational approximation given by VF, together with its deviation in terms of absolute error. Figure 14.Synthetic frequency response data and the rational approximation given by VF.

### 4.6. Single-phase transmission line modeling

Bode, Levy or LS, IRLS, VF, and LM are applied to the rational approximation of the frequency-dependent parameters corresponding to a single-phase transmission line. In Figure 15a, the diagram of the equivalent circuit for the test with L = 70 mH and Ri = 1.058 Ω is shown. The current source is considered as an ideal 60 Hz. The line length (l) is 100 km, with a diameter of 2.7 cm. The conductor is placed horizontally at a height of 20 m, with a 100 Ωm resistivity. The termination impedance at the end of line (Rl) is placed as 1e6 Ω, 1e-6 Ω, and 462 Ω for evaluation of the different line end cases. Figure 15.(a) Diagram of the equivalent circuit for the test, (b) behavior of Z(ω) and Y(ω).

The modeling consists of calculating rational approximations for the line series impedance Z and line shunt admittance Y by applying the aforementioned rational fitting techniques. Figure 15b shows the behavior of Z and Y as a function of frequency, calculated with 16,384 linearly spaced samples and Figure 16 shows a flowchart for the implementation of this test case.

Figure 17a shows the results for the fitting of Y with Np = 3, in terms of model deviations (relative error) by Bode, Levy, IRLS, VF, and LM. It can be observed that the best approximations are obtained with Levy, IRLS, LM and VF which converge to practically the same deviation. Figure 17.(a) Relative errors for the rational approximations of Y(ω) and (b) relative errors for the rational approximations of Z(ω).

Figure 17b shows the results for the fitting of Z with Np = 14, in terms of the model deviations (relative error) by Bode, Levy, IRLS, VF, and LM. It can be seen that the best approximation is obtained with VF, IRLS, and LM.

Then, the effect of the modeling errors in Y and Z on a time domain response is evaluated. The line terminal nodal admittance matrix Yn is established with respect to the two line ends 2 and 3 of Figure 15a as

I2I3=YnV2V3E70

where Yn is given by

Yn=YcAYcBYcBYcAE71

and Yc, A and B by

Ycω=Zω1ZωYωE72
Aω=cothZωYωlandBω=cschZωYωlE73

By using the Numerical inversion of Laplace Transform (NLT) [2, 12], the voltage responses on time domain (v1, v2 and v3) are calculated for the cases of open-ended, short-circuited, and perfectly matched lines end. These results are considered as a reference solution.

The rational function-based models obtained with Bode, Levy, IRLS, VF, and LM are used to calculate Z and Y. Using a similar procedure as mentioned for each technique, the voltage responses (v1, v2 and v3) are calculated through the NLT.

The objective is to calculate the error in the time domain introduced by the rational approximations of Z and Y in frequency domain, taking into account the reference solution.

The comparative results are shown in Figures 1820 for the cases of short-circuited, perfectly matched, and open-ended lines end, respectively. The absolute errors are consistent with the deviations for the rational approximation of Z. Figure 18.(a) Time domain simulation of the voltage response at the short-circuited line end, (b) absolute error with respect to the reference solution (v3). Figure 19.(a) Time domain simulation of the voltage response at the matched line end, (b) absolute error with respect to the reference solution (v3). Figure 20.(a) Time domain simulation of the voltage response at the open line end, (b) absolute error with respect to the reference solution (v3).

## 5. Conclusions

In this book, Bode, Levy, Noda, SK, IRLS, VF, LM, and DGN methods have been described in detail for complex-curve fitting, and a routine implemented in MATLAB environment presented for each technique. Rational approximation of synthetic frequency responses have used to show the operation of the programs. Moreover, the techniques are used to approximate the series line impedance (Z) and the shunt line admittance (Y) corresponding to a single-phase transmission line. The main conclusions are:

1. Asymptotic Approximation or Bode has proved to be a reliable technique for the modeling of overhead transmission lines. Also, this method can realize the fitting on the magnitude of the function and uses only real poles and zeros. However, the error level in the approximation can be substantial and it depends on the sensitivity criterion used when inserting new poles and zeros in the fitting.

2. The Levy method is pioneered for complex-curve fitting. Many techniques have been developed based on this methodology. Nevertheless, the method is biased, because it weights the fitting on high frequencies too much; this fact is demonstrated in the test cases.

3. The IRLS method is an accurate technique for the rational fitting of frequency responses. This method permits the implementation of different weighting functions in order to improve the level error in the approximation. Its disadvantage lies in the numerical ill-conditioning encountered in approximations with a wide frequency range.

4. The VF method is a robust and accurate technique. This fact has positioned this methodology as one of the most important in this field.

5. The LM method is a technique that shows good results in terms of level error. An advantage is that it can be implemented in pole-zero form, pole-residue form, and polynomial form.

6. The main disadvantage for the optimization techniques like LM and DGN is that they can get stuck in a local minimum.

7. The rational approximation for Z and Y in the single-phase transmission line modeling shows that the same technique does not always deliver the same result. Levy, LM and IRLS delivers more accurate result for the fitting of Y while VF reach the best result for the fitting of Z. The error levels obtained in time domain simulations are consistent with the fitting of Z, because in the modeling of overhead transmission lines, these parameters are more relevant, therefore VF reach the best level error.

compact PDF
Citations in RIS format
Citations in bibtex format

## More

© 2017 The Author(s). Licensee IntechOpen. This compact is distributed under the terms of the Creative Commons Attribution-NonCommercial 4.0 License, which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited.

## How to cite and reference

### Cite this compact Copy to clipboard

Eduardo Salvador Bañuelos-Cabral, José Alberto Gutiérrez-Robles and Bjørn Gustavsen (December 6th 2017). Rational Fitting Techniques for the Modeling of Electric Power Components and Systems Using MATLAB Environment, Rational Fitting Techniques for the Modeling of Electric Power Components and Systems Using MATLAB Environment, Eduardo Salvador Banuelos- Cabral, Jose Alberto Gutierrez-Robles and Bjorn Gustavsen, IntechOpen, DOI: 10.5772/intechopen.71358. Available from:

### Related Content

Next compact

#### Rational Fitting Techniques for the Modeling of Electric Power Components and Systems Using MATLAB Environment

By Eduardo Salvador Bañuelos-Cabral, José Alberto Gutiérrez-Robles and Bjørn Gustavsen

First chapter

#### Synthesis Processes for Li-Ion Battery Electrodes – From Solid State Reaction to Solvothermal Self-Assembly Methods

By Verónica Palomares and Teófilo Rojo

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.