The latest revised Fast-Relaxed Vector Fitting for MATLAB

This repo contains an updated and optimized version of the Fast-Relaxed Vector Fitting algorithm, originally developed by Bjørn Gustavsen. The present implementation maintains the conceptual fidelity of the original method, widely used in electrical engineering and dynamic systems analysis, but introduces some improvements in terms of computational efficiency, code structure, and compatibility with modern versions of MATLAB. The main objective of this version is to facilitate its use in current modeling and simulation projects by providing cleaner, more modular, and faster code. This version can be especially useful for researchers and practitioners who require a reliable, efficient and updated tool for analysis in the Laplace or frequency domain.

IMPORTANT NOTE 1: Download the VECFITX package (last release) and run the examples in the "example_cases" folder for testing. The vecfitX function is located in the "functions" folder, along with additional functions used internally by vecfitX.

IMPORTANT NOTE 2: If the equations in this README are not displaying correctly, please use Google Chrome.

Table of Contents

• 1. Vector Fitting Theory
  • 1.1. Classical Vector Fitting
  • 1.2. Improvement in pole relocation
  • 1.3. QR algorithm for efficient implementation
• 2. A comprehesive vecfitX tutorial
  • 2.1. New features
  • 2.2. Settings
  • 2.3. Function call
  • 2.4. Input and output data
    • 2.4.1. Description of the input data
    • 2.4.2. Description of the output data
  • 2.5. Iterative implementation
  • 2.6. Test cases
• 3. Contact info
• 4. References

Vector Fitting Theory

Classical Vector Fitting

One of the most utilized rational approximation techniques is Vector Fitting (VF). In 1999, VF was introduced and proved to be a highly robust and efficient rational approximation method, applicable to both smooth and resonant responses [1]. This caused VF to be rapidly adopted in several areas including power systems and macro-modeling systems. VF approximates a frequency response (generally an array) using a rational function approximation, as shown in (1):

where f(s) may represent the transfer function or the frequency response of a system, c_n and a_n represent the sets of residues and poles, respectively; and both d and e are real coefficients. Since the solution of (1) is a nonlinear problem, VF resolves this issue by performing the approximation in two stages, represented by two separated linear problems, as described next.

In the first stage, a set of predetermined poles â_n distributed over the frequency range of interest is used. If f(s) is a smooth function, real poles are employed; otherwise, complex poles are assumed, as in (2) :

Additionally, a frequency-dependent scaling function σ(s) is introduced with the same initial set of poles, as in (4):

Considering the multiplication between f(s) and σ(s), the following applies:

Substitution of (4) into (5) results in:

Since the initial set of poles is known, (6) represents a linear problem which can be solved via linear least squares approximation and provides the set of residues ĉ_n, which are used to obtain the zeros of σ(s) through the eigenvalues computation of (7):

where A is a diagonal matrix containing the set of starting poles â_n, b is a column vector of ones, and ĉ^T is a row vector containing the set of residues ĉ_n. Subsequently, (6) can be represented as (8) which leads to (9). Hence, the zeros of σ(s) are considered as the new poles of f(s).

where ẑ_n are the zeros of σ_fit(s), and z_n are the zeros of (σf)_fit(s).

In the second stage, VF uses the new set of calculated poles and (1) is solved in the least squares sense to obtain an updated set of residues ĉ_n. The aforementioned process can be performed iteratively until the error between the measured function and fitted function is minimum. The diagram summarizing this procedure is shown below:

Improvement in pole relocation: Relaxed condition of scaling function

VF operates by iteratively adjusting an initial set of poles to more optimal locations. When fitting the frequency-domain response of a rational function with the correct order, the poles can often be positioned accurately in a single step. However, when a lower-order function is used for the fit, multiple iterations may be required. The process becomes more challenging when the frequency response includes nonrational elements, such as noise, which can hinder convergence and, in some cases, even cause it to stall. A significant improvement in VF’s convergence properties can be achieved through a minor modification. As shown above, the classical formulation of VF incorporates a scaling function σ(s) that approaches unity at high frequencies. However, this high-frequency asymptotic condition can negatively impact convergence. To address this issue, the asymptotic condition is replaced with a more flexible constraint that ensures a nontrivial solution to the least-squares problem while avoiding the adverse effects on convergence [2]. Then equation (4) is replaced by (10):

where d̂ is real. In order to avoid the trivial solution. One equation is added to the resulting LS problem:

Equation (11) ensures that the sum of the real part of σ(s) over the given frequency samples is a nonzero value, while leaving all free variables unrestricted. Since σ(s) during iterations does not approach unity at high frequencies, (7) must be replaced with:

The zero calculation in (12) is only valid when d̂ is nonzero. If the absolute value of d̂ is found to be smaller than tol = 1x10^-8, the solution is discarded, and the least-squares problem has to be solved again with a fixed value for d̂ in equation (10).

QR algorithm for efficient implementation

For multiport systems, the classical VF method requires solving large, sparse linear systems to estimate the rational function, as shown in equation (6). This becomes computationally expensive and memory-intensive as the number of ports grows. For example, a 60-port system with 101 frequency samples requires 54 GB of RAM in the standard implementation, making it impractical for typical workstations. Even if the structure under study has a moderate amount of ports, the size of the corresponding LS matrix may become prohibitively large. Furthermore, a lot of computational effort is wasted on the calculation of the residues c_n, which are discarded by the VF algorithm. In [3], the authors introduce a Fast VF method that leverages the QR decomposition to simplify the LS equations. This process leads to a simplified set of equations which depend only on ĉ_n. Instead of solving one large linear system for all ports, the method processes each matrix element (port) sequentially. For each element, the QR decomposition is applied to its corresponding linear system. This decomposes the problem into smaller, manageable parts. After decomposing all individual systems, the results are combined into a smaller, shared system.

As explained above, once the residues ĉ_n of σ(s) are computed, the common poles â_n of the transfer function are found by computing the eigenvalues of (12). Note that the second stage of VF remains unchanged.

The modification of the VF with the two methods proposed in [2] and [3] have resulted in the third and last version of VF so far, better known as Fast-Relaxed Vector Fitting (FRVF) [4]. This version is implemented in this repository under the name vectfitX.

A comprehesive vecfitX tutorial

New features

The vectfitX routine introduces several improvements over the original implementation (vecfit3), making it more versatile and efficient for practical applications. Key features include:

• 3D Data array support: The function now accepts a 3D array as input for the matrix f(s). Unlike the original implementation, which required manual reshaping of the data into a 2D array (with frequency samples along the rows and matrix elements stacked in a single column), vectfitX automatically handles the input organization. Any 3D shape is accepted (e.g., rectangular, square, row or column vectors).
• Efficient handling of symmetric matrices: Symmetric data matrices are processed by fitting only the lower triangular part, significantly reducing computational effort. This was not supported in the original version, where users had to apply symmetry considerations manually.
• Cleaner and more readable code structure: The source code has been refactored to improve clarity, maintainability, and modularity.
• Improved computational efficiency: Redundant operations have been removed, and key processes have been optimized through compact expressions and vectorization.
• Future-proof MATLAB compatibility: The code has been reviewed and updated to ensure compatibility with current and upcoming MATLAB releases.
• Direct output of pole-residue representation: Unlike the original version, vectfitX directly returns the pole-residue model, eliminating the need for additional functions to extract this representation.
• Enhanced configuration options: Users can now customize output behavior more extensively. For example, generated plots can be automatically saved in common formats such as PDF, PNG, or JPEG.
• Improved plot aesthetics: The appearance of plots has been refined for better readability and visual presentation.

Settings

Like the original function written by Bjørn Gustavsen, vecfitX can be configured according to the user's requirements. Below is a list of the available options for customizing the function. You can set all the parameters or just some of them. Numbers in parentheses indicate the possible values ​​that a parameter can take.

• opt.relax:
  • (0) Use classical nontriviality constraint.
  • (1) Use relaxed nontriviality constraint.
• opt.stable:
  • (0) Unstable poles are kept unchanged.
  • (1) Unstable poles are made stable by flipping them into the left half-plane.
• opt.asymp:
  • (1) D and E are omitted in fitting process.
  • (2) Only E is omitted in fitting process.
  • (3) D and E are taken into account.
• opt.skip_pole:
  • (1) The pole identification part is skipped. C, D and E are identified using the initial poles as final poles.
• opt.skip_res:
  • (1) The residue identification part is skipped. Only the poles are identified.
• opt.repre:
  • (0) The returned state-space model has real elements only. Output variable A is square with 2x2 submatrices as diagonal elements.
  • (1) The returned state-space model has real and complex conjugate parameters. Output variable A is diagonal and sparse.
  • (2) The returned model has a residue-pole representation. Output variable A (poles) is a Nx1 vector, variable C (residues) is a Nrx(NxNc) array. Variables D and E are NrxNc arrays.
• opt.errplot:
  • (1) Include deviation in magnitude and phase angle plots.
• opt.fitplot:
  • (1) Create plots of fitted function compared to the original function. Both magnitude and phase angle are shown.
• opt.sigmaplot:
  • (1) Create plot of sigma function.
• opt.savefig:
  • (0) Figures are not saved.
  • (1) Save plots in PDF format.
  • (2) Save plots in PNG format.
  • (3) Save plots in JPEG format.
  • (4) Save plots in SVG format.

If you omit the vecfitX.m configuration, the function has the following default parameters:

def.relax = 1;                            % Use vector fitting with relaxed non-triviality constraint
def.stable = 1;                           % Enforce stable poles
def.asymp = 3;                            % Include both D and E  
def.skip_pole = 0;                        % Do not skip pole identification
def.skip_res = 0;                         % Do not skip identification of residues (C,D,E) 
def.repre = 1;                            % Create complex state space representation
def.errplot = 1;                          % Include deviation in magnitude and phase angle plot
def.fitplot = 1;                          % Create plots of fitted and original functions
def.sigmaplot = 0;                        % Exclude plot of sigma function
def.savefig = 0;                          % Figures are not saved

Function call

The following are some valid calls to the vecfitX function. The inputs and outputs are explained in detail below. It is important to note that vecfitX relies on several non-built-in functions. These functions, along with vecfitX, are included in the "source" folder. Please download the entire folder to ensure vecfitX functions properly. Test cases can be found in the "test_cases" folder. In some instances, external .txt files are required, and these are also included in the test cases folder.

[SER,poles,rmserr,fit] = vectfitX(f,s,poles,weight,opt);
[SER,poles,rmserr,fit] = vectfitX(f,s,poles,weight);
[SER,poles,rmserr] = vectfitX(f,s,poles,weight);

Input and output data

Description of the input data:

• f: Matrix f(s) function (3D array as shown in fig. 1.1) to be fitted with dimensions Nr X Nc X Ns. It is important to note that the elements of a layer (2D array) are listed by column (column-by-column order). vecfitX uses this order by default throughout the FRVF routine, even if a row vector is processed. In figure 1.1, f(s) is a non-symmetric matrix and is not limited to a square shape. In case f(s) is a symmetric matrix, vecfitX uses the lower triangular portion by default (as shown in Figure 1.2) to make fitting more efficient.
  
  

  Figure 1.1: f(s) representation as a 3D array (non-symmetric case).

  Figure 1.2: f(s) representation as a 3D array (symmetric case).

  • Nr: Number of rows in array.
  • Nc: Number of columns in array.
  • Ns: Number of layers (frequency samples) in array.
  
  
  Although f(s) is a 3D array, a 2D matrix representation is required within the FRVF algorithm. Before executing the algorithm, at the start of the vecfitX function, the original 3D array is converted into a 2D matrix using a column-by-column order (as shown above). In this way, the rows of the new matrix correspond to the number of elements in a layer (in the original array), and the columns correspond to the frequency samples of each element (as shown in fig. 2.1). In the case of a symmetric array, the number of rows in the new matrix is ​​lower because only the lower triangular part is used (see fig. 2.2).
  
  

  Figure 2.1: f(s) representation as a 2D array (non-symmetric case).

  Figure 2.2: f(s) representation as a 2D array (symmetric case).




• s: Vector of frequency samples [rad/sec] with dimensions 1 X Ns.



• poles: Vector of initial poles [rad/sec] with dimensions 1 X N. The selection of the initial poles can be done in different ways. Some examples of how to do this are shown in the test cases. Further details are given in [4].



• weight: The elements in the system matrix are weighted using this array. It can be used for achieving higher accuracy at desired frequency samples. If no weighting is desired, use unitary weights, i.e. weight array of ones, with dimensions 1 X Ns. Otherwise, 1D and 2D arrays are allowed.




Figure 3.1: Common weight array.

Figure 3.2: Individual weight array.

• 1D array: Common weighting for all elements, weight array with dimensions 1 X Ns. As shown in fig. 3.1, a multiplication is performed element-by-element between the weight array and each row of the reshaped f(s).
• 2D array: Individual weighting, weight array with dimensions (Nr*Nc) X Ns. In this case, the weight array has the same dimension as the reshaped f(s). Therefore, multiplication is done element by element between both arrays. If the original f(s) array is symmetric, the number of rows in weight array must have the same dimension as the matrix in fig 2.2. In both cases (symmetric and non-symmetric), before using any weighting strategy based on the elements of f(s), it is necessary to use stackM to transform f(s) into a 2D matrix as already explained. Different weighting strategies are shown in test cases. Refer to [4] for more details.



• opt: Configuration options.

Description of the output data:

• SER.A: Array containing the poles of the fitted matrix function.
  • If opt.repre = 0: A is sparse, square and real with 2 X 2 submatrices as diagonal elements, with dimensions N X N.
  
  
  

  Figure 4: Array A (SER.A) in the real state space model.

  
  
  
  • If opt.repre = 1: A is sparse, diagonal and complex, with dimensions N X N.
  
  
  

  Figure 5: Array A (SER.A) in the complex state space model.

  
  
  
  • If opt.repre = 2: A is a complex vector with dimensions N X 1.
  
  
  

  Figure 6: Array A (SER.A) in the residue-pole model.

  
  
  
• SER.B: Nx1 matrix only for real and complex state space representation.
  
  
  

  Figure 7.1: Array B (SER.B) in the real state space model.

  Figure 7.2: Array B (SER.B) in the complex state space model.

  • If opt.repre = 0: B is a column vector of zeros, ones and twos (see fig. 7.1).
  • If opt.repre = 1: B is a column vector of ones (see fig. 7.2).



• SER.C: Array containing the residues of the fitted matrix function.
  
  
  

  Figure 8.1: Array C (SER.C) in the state space model (real or complex).

  Figure 8.2: Array C (SER.C) in the residue-pole model.

  • If opt.repre = 0: C is real with dimensions (Nr*Nc) X N (see fig. 8.1).
  • If opt.repre = 1: C is complex with dimensions (Nr*Nc) X N (see fig. 8.1).
  • If opt.repre = 2: C is complex with dimensions Nr X (N*Nc) (see fig. 8.2).



• SER.D: If the residue-pole model is selected, D is a real constant term of dimensions Nr X Nc. Otherwise, D is a column vector of dimensions (Nr*Nc) X 1. Note that in a symmetric matrix the dimension of D is smaller, since only the lower triangular part is fitted.



• SER.E: If the residue-pole model is selected, E is a real proportional term of dimensions Nr X Nc. Otherwise, E is a column vector of dimensions (Nr*Nc) X 1. Note that in a symmetric matrix the dimension of E is smaller, since only the lower triangular part is fitted.



• ord_zrs: 1 X N matrix containing the newly calculated poles.



• rms: root-mean-square error (scalar) of approximation for f(s). 0 is returned if skip_res = 1.



• fit: (Nr*Nc) X Ns array containing the rational approximation of f(s). 0 is returned if skip_res = 1. If f(s) is a symmetric matrix, then fit has (Nr*(Nr+1)/2) rows and Ns columns. Note that fit has the same shape as the arrays in figures 2.1 and 2.2 (depending on the case). If you require a 3D representation of fit, just like the original f(s) function, use the res2fit function after calling vecfitX with the residue-pole representation, as shown below.

opt.repre = 2;
[SER,poles] = vectfitX(fs,s,poles,weight,opt);
pl = 1; % If pl = 1, res2fit creates plots of fitted function compared to the original f(s). Both magnitude and phase angle are shown.
[fs_fit, rms] = res2fit(s, SER, fs, Nr, Nc, pl); % fs_fit is the fitted function with the same 3D representation as f(s).

Iterative implementation

A common practice is to iteratively apply the FRVF to improve the relocation of the calculated poles and achieve higher accuracy. This is achieved as follows:

opt.skip_res = 1;          % Residue identification part is skipped.
opt.repre = 2;             % Pole-residue representation.
opt.fitplot = 1;           % Include graphs of fitted function.

Niter = 5;                 % Number of the iterations.
for iter = 1:Niter
  if iter == Niter
      opt.skip_res = 0;    % In the last iteration residue identification part is NOT skipped.
  end
    [SER,poles,rmserr,fit] = vectfitX(bigY,s,poles,weight,opt);
end

It is important to note that residue identification is skipped until the last iteration. This is to reduce computation time. Keep in mind that residue identification should be performed with the last set of poles calculated. When residue identification is disabled, graph creation is also disabled. Even if it has been manually enabled, as in the previous example. This way, only the plots corresponding to the last iteration are displayed.

Test cases

In order to demonstrate the use of the vecfitX function, test cases have been included in this repository. The examples and data have been taken and adapted from [4]. A brief description of these cases is shown below:

• ex1.m: In this program, an artificial scalar function f(s) is constructed. f(s) is a row vector whose elements correspond to two predefined partial fractions. In addition, 18 complex starting poles are used, vecfitX is called three times iteratively (new poles are used as starting poles in each iteration), and a real state-space representation is selected. The plots obtained are shown below.
  
  

  Figure 9.1: Magnitude plot.

  Figure 9.2: Phase angle plot.

  
  
  
• ex2.m: In this example, the measured response of a transformer f(s) (scalar) is read from file 03pk10.txt. In addition, 30 complex starting poles are used, vecfitX is called five times iteratively (new poles are used as starting poles in each iteration), a complex state-space representation is selected, and f(s) is weighted with the inverse of its magnitude. The plots obtained are shown below.
  
  

  Figure 10.1: Magnitude plot.

  Figure 10.2: Phase angle plot.

  
  
  
• ex3.m: In this program, a calculated terminal admittance 6 X 6 X 300 matrix Y(s) of a power system is read from file fdne.txt. Only the first column of Y(s) is fitted, as follows: 50 complex starting poles are used, vecfitX is called five times iteratively (new poles are used as starting poles in each iteration), a complex state-space representation is selected, and Y(s) is weighted with the inverse of the square root of its magnitude. The plots obtained are shown below.
  
  

  Figure 11.1: Magnitude plot.

  Figure 11.2: Phase angle plot.

  
  
  
• ex4.m: In this example, the same 6 X 6 X 300 matrix Y(s) is used. Now a complete complex state-space representation is calculated by fitting Y(s) column by column (each column has a different set of poles). In addition, 50 complex starting poles are used and improved on by fitting the element sum of the first column before fitting all columns, vecfitX is called 3 times iteratively (for each column), and the columns of Y(s) are weighted with the inverse of the square root of its magnitudes. The plots obtained of the last column are shown below.
  
  

  Figure 12.1: Magnitude plot.

  Figure 12.2: Phase angle plot.

  
  
  
• ex5.m: Again the same 6 X 6 X 300 matrix Y(s) is used. Now a residue-pole representation is calculated by fitting Y(s) (all elements has a common set of poles). The fitting is performed similarly to the previous example. After obtaining the residue-pole representation, the res2fit function is used to obtain the fitted matrix in its original shape. The plots obtained are shown below.
  
  

  Figure 13.1: Magnitude plot.

  Figure 13.2: Phase angle plot.

  
  
  

Contact info

If you would like to get in touch for professional collaboration, research proposals, or questions related to my work, feel free to reach out through any of the following channels:

• Primary email: josejamangape@gmail.com
• Secondary email: esteban.enriquez@cinvestav.com
• LinkedIn: Jose Esteban Enriquez
• GitHub: Esteban Enriquez
• ResearchGate: J. E. Enríquez-Jamangape

References

[1] B. Gustavsen and A. Semlyen, "Rational approximation of frequency domain responses by Vector Fitting", IEEE Trans. Power Delivery, vol. 14, no. 3, pp. 1052-1061, July 1999.

[2] B. Gustavsen, "Improving the pole relocating properties of vector fitting", IEEE Trans. Power Delivery, vol. 21, no. 3, pp. 1587-1592, July 2006.

[3] D. Deschrijver, M. Mrozowski, T. Dhaene, and D. De Zutter, “Macromodeling of Multiport Systems Using a Fast Implementation of the Vector Fitting Method”, IEEE Microwave and Wireless Components Letters, vol. 18, no. 6, pp. 383-385, June 2008.

[4] B. Gustavsen, "User's Guide for vectfit3.m (Fast, Relaxed Vector fitting)", SINTEF Energy Research, N-7465 Trondheim, Norway, 2008. Aviable online: vecfit3.m for MATLAB.