A Harmonic Balance Methodology for Circuits with Fractional and Nonlinear Elements

This paper discusses the ability to obtain periodic steady-state solutions for fractional nonlinear circuit problems. For a class of nonlinear problems with fractional derivatives (based on the Caputo or Riemann–Liouville definitions), a methodology is proposed to derive equations representing the dependencies between the harmonics of the sought variables. Two approaches are considered for how to address the apparent nonlinear dependencies: one based on symbolic computation and the other a numerical approach based on the analysis of time functions. An example problem with fractional and nonlinear elements is presented to illustrate the usefulness of the proposed methodology. Two error criteria are introduced to verify the accuracy of the obtained results. The methodology is mainly designed to provide referential solutions in analyses of the numerical method called SubIval (the subinterval-based method for computation of the fractional derivative in initial value problems).

Many theoretical considerations have been conducted concerning the behavior of fractional-order systems (e.g., stability and controllability analyses [27]). The most important theoretical aspect is, however, the ability to solve problems where fractional derivatives and integrals appear, as this is essential for all other analyses.
Various definitions of the fractional derivative can be found in the literature [25], but the most commonly used definitions are those of Riemann-Liouville [36] and Caputo [8].
This paper considers circuit analyses, specifically those with periodic steady-state sources and fractional or nonlinear elements, even elements that are both fractional and nonlinear. The motivation for the research is given in the next section.

Motivation
The current paper presents part of the studies concerning the design of numerical and semi-analytical methods (providing referential solutions) to solve problems with fractional derivatives. The study mainly concerns problems emerging in circuit analyses.
When designing a numerical method, it is appropriate to have means for a reliable assessment of its usefulness. The efficiency of a method is often estimated through analyses of errors [30,58] and computational components (e.g., the number of basic operations of an algorithm [54] or computation time [33]). When a method has already been implemented, it is worthwhile to test it on selected problems, especially those for which the method is mainly intended. The correctness of the solutions obtained through the method can be determined through: -fulfillment of equations specific for the problem, e.g., basic laws, such as Kirchhoff's laws or the power balance for a circuit problem [39]; -comparison with results obtained through another method, preferably one operating on a very different basis (e.g., the results obtained through the application of numerical methods are compared with results from analytical solutions [35] and vice versa [53]).
The current study was motivated by the author's work on the numerical method called SubIval [52] (the subinterval-based method for fractional derivatives in initial value problems). The method is designed mainly with circuit problems in mind, but it is not limited to them in its application.
While determining the accuracy of the method for solving circuit problems, it has been established that: -for selected transient problems, the solutions can be compared with results obtained through the evaluation of analytical solutions based on the Mittag-Leffler function [18,22,23,34]; -steady-state linear AC problems solved with the newly designed numerical methods can be compared with solutions obtained through the application of complex numbers.
It would also be useful to have the ability to obtain referential solutions for nonlinear problems. As suggested previously, it would be ideal if the method for obtaining these solutions had a basis that is very different from that of the method being examined.
The presented methodology has been designed for periodic steady-state nonlinear problems of circuit theory, where the time functions (of which the solution comprises) consist of limited numbers of significant harmonics. The methodology leads a problem to a harmonic balance form. Optionally, the methodology can be based on symbolic computation when acquiring dependencies that result from nonlinearities. The main goal is to obtain steady-state solutions for a selected class of problems that can be used for further improvement in the SubIval numerical method. The considered class of problems is one that results in Eq. (8) described in Sect. 4.

Assumptions and Complex Number Representation
The current study is simplified to the case where: Assumption (c) leads to an obvious source of error as harmonics above h max are removed from the result, whereas in reality, even higher harmonics will emerge. This motivates the "harmonic remainder error" defined in Sect. 9.
Every time function of the solution (for the moment denoted by w(t)) can be given by a vector of complex numbers representing the subsequent odd harmonics: Hence: where the parameter: with f 1 being the base frequency.
The study considers the application of the Riemann-Liouville [36] and Caputo [8] definitions of the fractional derivative, where the order is α ∈ (0, 1]. For a steady-state analysis, both definitions result in the same operation because they differ by only a component dependent on the initial value [2]. The fractional derivative (in the periodic steady state) can therefore be computed by: where: For a fractional integrodifferentiation in a complex vector representation (resulting in a vector d), one can use the Hadamard product: where:

Time-Dependent Form
The methodology is designed for problems that yield a system of equations in the following form: where x(t) is a vector of n x state variables, and y(t) is a vector of all the remaining n y variables. The vector v(t) (of size n v ) contains all the source time functions. The vector w(t) contains all the variables being computed: hence, the total number of computed variables n w = n y + n x . The component D α t x(t) contains the fractional derivatives (each in either the Riemann-Liouville definition or Caputo definition): with α denoting the vector of fractional derivative orders. As for the other components: -M I is an n y × n y matrix; -M II is an n y × n x matrix; -M III is an n x × n y matrix; -M IV is an n x × n x matrix; -T is an n y × n v matrix; -F NL (w(t)) contains n NL nonlinear functions (generally dependent on w(t), but in truth each depends on only one variable); -a vector denoted by 0 k is one consisting of k zeros.
Note that nonlinearities are given in the form: where i LHS is the index of the left-hand side variable, which depends on f NL i . i arg is the index of the argument of the nonlinear function. Hence, the left-hand sides contribute to entries in M I or M II , and each nonlinear dependency is stored in F NL (w(t)).
For further convenience, the nonlinear dependencies are represented as follows: where it is assumed that the dependencies on the appropriate variables are indexed by integers placed in the auxiliary vector i arg (the notation (i arg ) k denotes the kth element of the vector i arg ). Additionally, another auxiliary vector, called i LHS , has been used to store the indices of the left-hand side variables of the nonlinear functions in the form of (11). This application of auxiliary integer vectors has greatly aided the implementation of the methodology.

Conversion to the Harmonic Balance Form
The current subsection discusses the core of the proposed methodology. For a periodic steady-state problem resulting in a system of nonlinear equations in the form of (8) (including fractional differential equations), the proposed method provides the ability to obtain an alternative system of nonlinear equations representing relations of the sought variables' harmonics. The solution is sought in the form of real and imaginary parts (i.e., the sine and cosine components of the actual time function) for each harmonic h of variable w i . This approach is common [11,29,55], although the magnitude and angle of the complex number could also be considered (as in [31]). The harmonics are denoted by The solution vector is: The system of nonlinear equations has the form: Subject to (13), matrix A can be rewritten in the following convenient form: where the auxiliary matrix: contains the multipliers of both the sine and cosine components of harmonic h of the solution [which is why each A h appears twice in (15)]. The component S h R results from (7); it is a diagonal matrix of the form: The matrix denoted by A j features only components that arise from (7); it can be presented as: where all A j h are sparse matrices of the form: where the notation 0 k× j represents a zero matrix with k rows and j columns. S h I are diagonal matrices that (like S h R ) result from (7). They have the following form: The right-hand side vector b in (14) can be filled after determining the harmonics of N(w) is a sparse vector filled only for entries resulting from the nonlinear equations of (8). Since each nonlinear function of F NL (w(t)) in (8) actually depends on only one variable of w(t), the dependencies can be written as: where the nonlinear dependency N i,h,part (w (i arg ) i ) is the appropriate part (sine or cosine) of the hth harmonic, resulting from the original nonlinear dependency f NL i (w (i arg ) i ). Two methods of addressing the nonlinear dependencies are discussed in Sect. 6.

Remarks on Nonlinear Solver
Because of the composite form of the resulting nonlinear dependencies given in N(w), which could introduce problems when attempting to solve the resulting nonlinear system of equations, a stage-based methodology is proposed. Such an approach, whether used for solving systems of equations or optimization problems [32], can be very useful when a large number of unknowns is considered and when there is uncertainty about what problems the nonlinearities can introduce for a selected solution technique.
In the first stage, a single-harmonic solution is assumed, and the resulting system of equations is solved. This solution is used as the starting point for the first time harmonic when seeking the solution in the next stage, i.e., for h max = 3. In each subsequent stage (up to the actual h max ), the solution of the previous stage is used as the initial guess for the harmonics up to h max − 2.
Furthermore, only the final stage requires a low error tolerance for the internal iterative solver being used. The error tolerance for all the other stages can be much smaller as the role of these stages is only to get close to the solution.

Handling Nonlinear Dependencies
Two approaches to addressing the nonlinear dependencies in N(w) are described in this paper. The first is based on symbolic computation, and the second relies on the extraction of harmonics from periodic time functions.

Symbolic Computation Approach
In this approach, it is assumed that the nonlinear dependencies are, or can be, approximated (with reasonable accuracy) by an odd power series with maximum power k max .
When assuming a complex harmonic representation, whenever a nonlinear dependency appears, it is first converted into a set of dependencies for the appropriate harmonics. A nonlinear function f NL is represented by an odd power series following assumption (b) of Sect. 3: The result has a maximum time harmonic h = k max h max . The result of each exponentiation of the time function w by the integer k, cut down to harmonics h = h max and below, can be represented by a collection of nonlinear dependencies: with w following Eq. (1).
To derive the dependencies in symbolic form, the author has applied a simple lightweight library for the symbolic computation of multivariate polynomials [48]. An efficient algorithm for multivariate polynomial multiplication greatly decreases the time needed to perform the required symbolic computations, even when considering many terms [50].
The base operation for obtaining the harmonics' dependencies is multiplication between two harmonics h = i and j ≤ i of the harmonics vectors w and u. This multiplication contributes to the harmonics h = i − j and h = i + j of the total result, as depicted in Fig. 1.
The real and imaginary parts of the nonlinear dependencies given by the symbol Υ generally have the form of multivariate polynomials in expanded form: where a i are real-valued multipliers, k represents nonnegative exponentiations, and N is the number of monomials.
As h max and k max increase, the dependencies take the form of very large symbolic objects, where not only memory but also the computation time could be an issue. However, one advantage is that these computations need to be performed only once for all integers k ≤ k max and a selected h max , and the results can be saved for later use. The formulation of any nonlinear dependency collection: (27) and other operations that could be required during the solution process (such as differentiation during Jacobian matrix evaluations) do not introduce much numerical weight for efficient symbolic computation implementations [49].
The dependencies in each f NL (w) can be used when addressing N(w) in (23).

Time Function Approach
The second proposed approach to address the nonlinear dependencies is based on a much simpler concept. Instead of knowing the form of N(w), as in the symbolic approach, the actual dependencies are handled only during evaluations of a nonlinear system solver.
The evaluation of N(w) is performed according to the following instructions: (c) for each of the points on the time axis, the time function of the nonlinear function's argument υ is evaluated, resulting in values denoted by υ 1 , υ 2 , . . . υ n t ; (d) the results (denoted by r j ) are computed: (e) the harmonics for h = 1, 3, . . . h max are obtained with a selected algorithm (e.g., fast Fourier transform or the Vaníček method [57]); these harmonics are then used to fill the vector resulting from the evaluation of N(w).
A comparison of some of the properties of the two approaches (with respect to how the nonlinear dependencies are handled) is given in Table 1.

Remarks on the Implementation
The conversion from (8) to (14), according to the harmonic balance methodology, is performed with programs and libraries written by the author in C#. The Math.NET Numerics library [42] is applied for matrix and vector operations.
The harmonic balance methodology can be applied for a problem in the form of 8. All the necessary information about the problem is put into an object, which is named prob. The base of the A matrix, which is independent of frequency (formulated from the M I , M II , M III and M IV matrices), can be formulated using the following fragment of code. The matrices are all Matrix<double> objects of the Math.NET Numerics library, and the vectors appear as Vector<double> objects. Next, the A + A j matrix and the b vector of (14) are prepared. Additionally, the nonlinear dependencies of the problem are copied and are then applied to enable evaluation of N(w) on the right-hand side of (14). The formulations must be done for all stages of the nonlinear solver described in Sect. 7. The following code fragment is applied. The NLinfo object contains information about the nonlinear dependencies of the original system. A separate procedure controls the evaluations of N(w), either using evaluations of symbolic expressions or following the time function approach described in Sect. 6.2.

Example 1
The discussed methodology has been verified via computational examples. Two of these examples are presented in this paper. The first concerns the circuit presented in Fig. 2. The circuit is purely theoretical, but it features elements that appear in analyses of real circuits, e.g., strictly increasing nonlinearities, which ensure a unique solution.
The fractional capacitor is described by the equation: where u cmn denotes the common voltage of the fractional capacitor, nonlinear coil and nonlinear resistor. For the fractional, nonlinear coil, the differential equation is: where ψ is the pseudo-flux of the fractional coil (unit: Wb · s γ −1 ). The nonlinear function ψ(i ψ ) is given by: where ψ 0 and i 0 are the nonlinear function parameters: The nonlinear resistor R NL is described by the equation: with the nonlinear function parameters: For the considered circuit, it is possible to formulate equations in the form of (8) and, by applying the harmonic balance methodology, to transform the problem to the form of (14).
The transformation of presented problem to the form of (8) is described in Appendix A.
For the two approaches to addressing the nonlinear dependencies, if the symbolic approach is selected, then the arctangent function must first be approximated by a power series. First, however, it is worthwhile to modify the dependency from ψ(i ψ ) to i ψ (ψ): because fewer power series terms are required. However, for the analysis presented in this paper, the time function approach is applied to directly address the original nonlinear dependencies.
The solution of the time periodic problem is obtained for h max = 25 in the form of the vector w, as given in Eq. (13). Then, the solution is converted to a time-dependent form. The results (for 3 periods) are depicted in Fig. 3 along with a comparison with a numerical solution obtained by means of a solver using SubIval [47,52].
The computations of SubIval use the author's C# programs and apply a DLL available at [47]. The library uses part of the code given in [10] to compute the gamma function. The SubIval step size adaptive solver requires a system of nonlinear equations to be solved at each iteration. For this purpose, a modified Gauss-Newton method, which applies methods and classes from the Math.NET Numerics library, is implemented.
The solver can obtain a time-dependent solution; hence, n T = 10 periods of the solution are obtained. Then, only the final solution is compared with the periodic steady-state solution obtained by solving the system of equations resulting from the harmonic balance methodology. All the initial conditions for the state variables (i L , u cmn and ψ) are set to zero, and the variables selected for comparison are the state variables.
One can notice that as the number of periods increases, the results become closer to each other. To reliably ascertain how close the results are, two criteria have been proposed, both of which are explained in Sect. 9.
The time required to complete this task 1 for h max = 5 is only 0.26 s, whereas h max = 25 requires 14.73 s. Naturally, the more complicated the problem is (especially when adding nonlinear elements), the more the computation time will increase when taking into account more harmonics, especially if the nonlinear system is to be solved with high accuracy (note that the final stage of the nonlinear solver is always executed with a lower tolerance for the objective function). Naturally, the methodology can also be applied in an integer-order case where β and γ are equal to 1. The manner in which the general equations for the problem [i.e., in the form of (8)] are formulated is the same, but when applying the harmonic balance methodology, the matrix denoted by S h R in Sect. 4.2 is empty and: because all the derivative orders are equal to 1. The SubIval solver can still be applied because SubIval itself supports first-order time derivatives [52]. The results of the comparison in the integer-order case are presented in Fig. 4.
In the case of integer-order derivatives, the waveforms also indicate a steady-state result similar to the numerical solution. The results are also verified through the criteria discussed in Sect. 9.

Example 2
The second example concerns the steady-state solution of the circuit presented in Fig. 5. This example is also purely theoretical.
The circuit features two nonlinear coils, where the first (denoted by the symbol L ψ ) is described by the same nonlinear function as in the previous example [this time denoted by ψ = ψ(i ψ ), given in Eq. (31)] and the differential equation: The second coil is described by the fractional differential equation: Φ is an artificial variable (referred to as the "pseudo-flux", whose unit is Wb · s λ−1 ) of the second coil. The relation between Φ and the current i Φ is described by the nonlinear equation: The parameters taken for the computation in this paper are: The circuit features one more element that is both fractional and nonlinear, the capacitor denoted by C q . The relation between its voltage u q and the variable q (with the unit C · s β−1 ) is described by the equation: where the following values are assumed: The differential equation for this element is: The nonlinear resistor R NL follows the same nonlinear dependency as in the example in Sect. 8.1 (in this case, denoted by i ψ (u NL )).
The formulation of a system of equations in the form (8) for the discussed problem is presented in Appendix B.
The solution for the problem is obtained for h max = 25, which was sufficient to obtain a solution with satisfactory accuracy (with respect to the criteria defined in Sect. 9). The solutions are again compared with those obtained through the time stepping solver applying SubIval. The final (n T th) period of the solution is selected for the error calculations. The state variable time functions (ψ, Φ, i L , u C and q) are selected for the comparison. The first three periods of the solutions are depicted in Fig. 6.
Because this problem contains more variables and more nonlinear dependencies than in the previous example, the computation times are longer. For h max = 5, the task required only 1.32 s; however, for h max = 25, the duration was much longer, i.e., 264.66 s.
For this example, the solution is also been verified through the criteria explained in Sect. 9.
As in the previous example, the harmonic balance methodology can also be applied for an integer-order case, where α is a vector of ones. The results for this case are presented in Fig. 7. The first criterion by which the accuracy of the solution is determined is comparison with results obtained through the SubIval solver, which can obtain highly accurate solutions for linear transient problems and AC problems [46,51,52]. The most important parameters for the SubIval solver are given in Table 2.
The above parameters are used for both problems (note that T is dependent on the base frequency, which is 50 Hz for the first problem and 60 Hz for the second problem).  Maximum of e an versus num 5 6 .91 × 10 −1 1.11 × 10 1 1.90 25 9.76 × 10 −3 1.59 × 10 −1 1.14 × 10 −1 p mov is the order of the polynomial approximation used in the core computations of SubIval, e max is the maximum allowed estimated error, e ctrl is the desired value of the estimated error (according to which the adaptive solver modifies the time step), and t min and t max are, respectively, the minimum and maximum values of the step size.
The results obtained for the harmonic balance methodology are compared with the numerical solution for time nodes t 1 , t 2 , . . . t n of the last period (i.e., t ∈ [(n T − 1)T, n T T ]). The time nodes are those selected by the SubIval solver during the time stepping process. The error for a selected variable w is computed according to the following formula: where w num i is the value obtained by the SubIval solver for the time instance t i .

Results for Example 1
The maximum and average of the obtained error values, for the selected variables (i L , u cmn and ψ), are given in Table 3. The error has been checked for the case of h max = 25 (mentioned in the previous section), along with a less demanding solution obtained for h max = 5.
The error values show that as the number of considered harmonics increases, the results get closer to those obtained by means of the SubIval solver. Table 4 shows the average and maximum error values for the case where the inertial elements are replaced by their integer-order alternatives.
Again, smaller error values are obtained for larger h max , which indicates that the methodology also works for the integer-order case and can be applied to classic periodic steady-state nonlinear problems in circuit theory. Additionally, the errors are in a similar range as those for the fractional case.   Table 5 presents the error values obtained for the fractional case of the second example, and Table 6 gives the results for the integer-order case.
As for the error values in the first example, more accurate solutions are obtained when h max is larger according to the considered criterion. Additionally, the same conclusion can be drawn for the integer-order case.

Harmonic Remainder Error
Clearly, for the h max harmonics, each nonlinear dependency will yield harmonics greater than this number (which are later not taken into account). This practice of cutting off harmonics is necessary when analyzing the result for each variable as coefficients of a Fourier series.
After the solution is obtained, one can estimate the harmonics that have been cut off for each nonlinear dependency, which can be done by evaluating the harmonics through multiplication (such as when the symbolic computation approach is selected to address the nonlinear dependencies). However, this process would require each of the nonlinear dependencies to be given as a power series. To directly address the nonlinear dependencies, the time function approach is selected for the computations described in this paper. The harmonics that have been cut off are determined in a similar way. This strategy can be viewed as simply the application of the time function approach with an attempt to obtain additional harmonics (h remainder > h max ) in the final step (step (e) described in Sect. 6.2). h remainder = 51 is selected for further computations. Through this approach, one can estimate not only the omitted harmonics but also a measure of correctness for the harmonics appearing in the obtained solution. Both the included and omitted harmonics in the solution are taken into account in what is later called the "harmonic remainder error".
For each nonlinear function, the harmonics of the right-hand side (denoted by R h each) are obtained for the mentioned harmonics up to h remainder . Assuming that L h denotes the harmonics of the variable on the left-hand side of the nonlinear equation [note that the nonlinear functions still follow the form of Eq. (11)], an error value can be computed for each of the harmonics: The harmonic remainder error is computed by the formula: where L max denotes the maximum value of the time function obtained for the left-hand side variable of the nonlinear equation (i.e., the maximum taken from an arbitrarily selected sufficient number of points).

Results for Example 1
The harmonic remainder error values obtained for h max = 25 and, once again, the less demanding h max = 5, for the nonlinear functions i NL (u cmn ) (of the nonlinear resistor R NL ) and ψ(i γ ) (of the fractional, nonlinear coil L ψ ) are given in Table 7. The harmonic remainder error values obtained for the integer-order case are given in Table 8.

Results for Example 2
The harmonic remainder error is been computed for the second example (presented in Sect. 8.2). The results of these computations (for the nonlinear functions ψ(i ψ ), i Φ (Φ), i ψ (u NL ) and u q (q)) are presented in Table 9. Again, the error is computed for two cases, i.e., h max = 5 and the more demanding case of h max = 25. The error is also computed for the integer-order time derivative case. The results (also for h max = 5 and h max = 25) are given in Table 10.

Summary
For the class of nonlinear problems (with fractional derivatives) resulting in the system of Eq. (8), a methodology has been presented, which for a sought periodic steady-state solution, allows the problem to be converted to a harmonic balance form [generally described by Eq. (14)]. In the new form, the unknowns are the coefficients of the sought variables' harmonics up to a selected h = h max . Two approaches have been proposed to address the nonlinear dependencies in the harmonic balance form. The first is a symbolic approach, whose advantage is that the computational weight is placed on a pre-solution process, the results of which can then be used conveniently (with a gain in efficiency when the solution is sought). The disadvantage of the symbolic approach is the requirement for each of the nonlinear dependencies to be given as a power series. The second approach (based on time functions) allows the original nonlinear dependencies to be used. Its disadvantage is that the entire computational weight is moved to the solution process.
Two computational examples have been presented to demonstrate the usefulness of the harmonic balance methodology. The steady-state solutions have been compared with results obtained through a selected numerical method-an adaptive time step size solver applying SubIval [46,51,52].
To verify the results, two error criteria have been introduced. The first one is based on comparison with the numerical result, and the second determines the accuracy by estimating the significance of the harmonics not included in the result.
In future papers, more examples will be presented. The harmonic balance methodology has proven to be efficient for solving fractional nonlinear circuit problems in periodic steady states. It has been mainly designed to provide reference solutions for the analysis and improvement in SubIval and solvers applying it.
The methodology has the potential to be very accurate because it allows solutions to be obtained without errors emerging from the estimation of the fractional derivative.
In future analyses, it would be interesting to determine whether the methodology can be applied with equal success to other definitions of the fractional derivative that are commonly applied in circuit analyses, such as the Atangana-Baleanu definition [4,20] and others [3,5].
This methodology is valid only for steady-state analyses. An idea for its extension could be to seek a more general solution that applies generalized cos α and sin α functions, as presented in [41]. A future study on such a method could fill the spot in Table 11 where non-time stepping methods are given, which can provide reference solutions when assessing the accuracies of numerical solvers for circuit problems.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A System of Equations for Example 1
This section explains how a system of equations in the form of (8) can be formulated for the first computational example presented in this paper (depicted in Fig. 8 along with auxiliary variables).
The vector of state variables is: with the derivative orders: while the vector of additional variables is: The voltage source time function constitutes the only entry in the source vector: From the dependency between the voltage source and the node potential V 1 , one obtains the entries: Fig. 8 The first example: circuit with integer-order coil, fractional-order capacitor (of order α), fractional (nonlinear) coil (of order γ ) and a nonlinear resistor. Additional auxiliary variables have also been marked the auxiliary integer entry: and the nonlinear dependency The remaining entries of (8) are those of M III and M IV , which determine the form of the differential equations. For the integer-order coil: For the fractional capacitor: Finally, for the fractional, nonlinear coil:

B System of Equations for Example 2
This appendix provides information on how a system of equations in the form of (8) can be formulated for the second example discussed in this paper. The example is presented in Fig. 9. The equations for the problem can vary depending on the selected auxiliary variables. The vector of state variables is: Fig. 9 The second considered example with fractional and nonlinear elements. Additional auxiliary variables have also been marked while the respective derivative orders are given by: while for the middle branch: