Optimization Techniques for Permittivity and Permeability Determination

This paper discusses optimization techniques for the determination of complex permittivity and permeability in transmission lines. The traditional theoretical model using scattering parameters is extended into a mathematical regression model that can be solved with widely accepted numerical techniques. This new model produces accurate primary mode results for the samples tested including nonmagnetic and magnetic materials with high dielectric constants. An extension of the model includes responses due to higher order modes. The general model determines parameters to specify the spectral functional form of complex permittivity and permeability and is capable of small corrections to independent variable data including angular frequency, sample length, sample position, and cutoff wavelength. The method provides reliable determination for both low and high permittivity materials.


Introduction
A constrained nonlinear optimization procedure is presented for the determination of complex permittivity and permeability spectra from scattering parameter (5-parameter) data taken from an automatic network analyzer (ANA). The procedure has been used successfully for reliable characterization of permittivity and permeability of many different test samples. In addition, it provides a basis for the analysis of multi-mode field data and for the determination of experimental systematic uncertainty.
Previous work in this area involved the determination of permittivity and permeability on a pointby-point basis with explicit or implicit solution of a system of nonlinear scattering equations at each particular frequency (see [1,5]). Inaccurate results, however, may arise when numerical singularities occur at frequencies corresponding to integer multiples of one half wavelength of the material. The same system of nonlinear scattering equations is used in this study. Here they are solved in the sense of least squares over the entire range of measurement by determining the best Laurent series approximations to permittivity and permeability consistent with linearity and causality constraints. Points of singularity may be de-emphasized to lessen the effect of highly uncertain data points.
This effort determines complex permittivity and permeability from two-port 5-parameter data using the primary (or fundamental) mode field behavior in various materials. Physical measurements of the 5-parameter data are made with an automatic network analyzer. Fundamental mode 5-parameter relationships are then used to solve for permeability and permittivity, /i.(<u) and e(w) respectively, as a function of the angular frequency <u.
The general functional form for ^ and e was decided upon after evaluating several different polynomial and trigonometric relations. The evaluation criteria for the functional form for permittivity and permeability was based on the reproducibility of the 5-parameter data in terms of the minimum total least square approximation to these data. The best overall approximation was selected. The general functional form for /x and e with acceptable results involves only the first two terms of a Laurent series.
The optimization approach in this research is an implicit function regression model. This model is solved by the orthogonal distance regression package ODRPACK [2]. This estimation package allows for adjustments in input parameters to compensate for measurement uncertainties. Adjustments here are limited to the sample length, sample position in the waveguide measurement fixture, and the cutoff wavelength. Orthogonal distance regression is intended to compensate for slight uncertainties in the independent variable (angular frequency, co) as well as the dependent variable (observed 5-parameter data).
The physical model is outlined in Sec. 2, with the related mathematical model discussed in Sec. 3. Numerical considerations are covered in Sec. 4 and conclusions and future directions discussed in Sec. 5.

Scattering Parameter Relations
The equations described below relate the measured two-port scattering parameters (5-parameters) to the permittivity and permeability of the material. First, in order to develop the scattering equations, the following notation is used. Let be the corresponding angular frequency. Then the permittivity and permeability of a sample material, where eo and JLIO are the permittivity and permeability of a vacuum, and ER and IIR are the relative complex permittivity and permeability.
Next let Cvac and ciab be the speed of light in a vacuum and the laboratory, respectively, and for a given frequency/, let It is assumed that the total length of the sample holder is (1) as shown in Fig. 1, where Li and Li are the distances from the calibration reference planes to the sample faces for ports 1 and 2, respectively.
For a two-port device the expressions for the measured 5-parameters are obtained by the solution of a related boundary value problem. The explicit expressions for the scattering relations of the fundamental mode are assumed to functions of \ci are given by 52i(Ac,) =Ri Ri I i-fi^ij' where Ri = exp( -yo Li) and Ri = exp( -yo L2) are the reference plane transformation expressions.

IJ-K(CII), e*(io), and the Regression Model
The mathematical problem can be stated as that of finding parameters to a prespecified functional form for the complex functions /Xfi(w) and e^(w). S-parameter data acquired from the ANA for selected samples of materials provide complex values for each of Sn, S21, Sn, and 522, at n different frequencies ranging from 1 to 18 GHz. The model determines the parameterization of JXR and e^ that best reproduces simultaneously the four 5-parameters for the n observations, given the sample data, two reference plane positions, and the sample length.
The general form for /i«((o) and €«(«) uses terms from the Laurent series: /(«)= 1 Oi (1-1-6,0))'' (6) where a, and bi are complex scalars in the /th term. The solution procedure here is an implicit function regression and corresponds to the parameterization of two terms in a truncated Laurent series. Note that Eq. (7) automatically satisfies Kramers-Kronig relations (see, for example, [1] or [4]) for dispersion. The functional form for ^.^(w) and e^(w) follows: For notational convenience, the terms of the second truncated Laurent series ^1,^2, Bi, and B2 will be referred to as Az, A4, B3, and B4.
The solution procedure to determine the complex parameters At, Bi is the minimization of the sum of the squares of the uncertainties between the predicted and observed 5-parameters, where, for ij = 1, 2, 5,* represents the fcth observed Sij scattering parameter at frequency w* and

Adjustments to the Model Inputs
In standard ordinary least squares regression models, the observed responses are assumed to contain some uncertainties either produced by the phenomena under examination or introduced by the device that measures the events. In addition, certain independent and dependent variable pairs may not be as reliable as others due to an increased variance in the uncertainties in the dependent variable for particular values of the independent variable. The application of a reduced weight to points of questionable reliability provides the modeler with a means to de-emphasize such points to find a more appropriate regression solution.
Usually the independent variable can be controlled, and the precise value of each of the observations is well known. An orthogonal distance regression model provides the modeler with the additional ability to assume that the independent variable, in this case frequency, may contain some uncertainty as well. Allowances for this type of uncertainty can, in some cases, greatly improve the approximation. For this particular model and the samples tested in this study, the uncertainty in the independent variables is sufficiently small to allow the modeler to assume that an ordinary least squares approximation provides an adequate solution.
Other model parameters such as sample length, sample position in the waveguide, and cutoff wavelength are sufficiently sensitive to require slight perturbation. Each of these inputs is required by the 5-parameter equations and is considered to be known by the modeler. Since these inputs are not always known exactly, each may be perturbed slightly, as determined by the problem solver, to improve approximation. This allows the user to adjust for measurement uncertainty.
In particular, since incorrect specification of the sample position, Li, in waveguide affects the value of phase in the reflected 5-parameter data, an additional parameter, /3i.i, is included inRi, Similarly, uncertainties in Lair and L2 are represented by analogous means with parameters jSm, and /3i2, respectively. With Eq. (1), the total length L of the sample is completely determined by and is parameterized by the values of PL^" PLI, and Another parameterized correction, fix, is included in the cutoff wavelength A^, as Ac i + JSA . Inaccuracies in the waveguide dimensions due to the milling process can affect the value for the cutoff wavelength. In addition, the higher order model discussed below requires an additional wavelength cutoff value, X^^, for a higher mode solution. Since it is not known a priori which higher modes will be present, a variation in \c2 from 0<Xc2^ Ki + I^x, enables the solution procedure to find the value for Ac2 that best improves the 5-parameter approximation for the higher order terms.
One can formulate an invariant model with respect to reference planes for the problem discussed here, as suggested in [1], to remove uncertainty in L1 and in L2. This approach uses Eq. (3) or (4), and the determinant of the 5-matrix, det(5)=5ii522 -5i252i, and solves for the parameters Ai,..., AA, Bi ..., 54. Simplification of det(5) yields a formula with neither Li nor Li. This reduced model remains dependent on a precise knowledge of both Lair and L and the cutoff wavelength, Ac,. This approach can be used interchangeably with the original model that contains Li and L2. The original model uses twice the number of observations and seems to produce more accurate approximations to the 5parameter data.

Higher Order Modes
In samples with a high dielectric constant, the observed 5-parameter data may exhibit responses due to modes other than the fundamental mode. These responses are the result of resonances of the higher mode in the material. Earlier work does not include this information in the computation of /n^ and ER as higher modes are ignored. In this section we describe an enhancement to the previously defined model that does include higher order response data.
Because of the similarity between the response of the primary mode and the higher modes (see Sec. 4.4), an additional term is added to each of the four equations (see Eqs. [2][3][4][5]. This numerical model includes higher order mode structure in the 5-parameter approximation and is a simple extension of the primary mode model described in Sec. 1} The terms that approximate the higher order modes are identical to the primary mode term with the exception that each term is scaled. For the evaluation of the higher order terms, all input parameters are unchanged except for the parameter for the higher mode cutoff wavelength, Ac 2, that is allowed to decrease. The explicit forms of the new set of equations use the Eqs. (2)(3)(4)(5), and are denoted as Dw, Dn, D21, and D22. Associated with each higher order term is a scaling parameter, )8i, ..., pn. In particular, to include the primary mode and one higher-order mode, the following approximate model is used: I>12=5l2(AcO + ^35l2(AcO, (11) £>22 = 522(Ac.) + p4S22{\c2), where Ac^ is determined by ODRPACK.

The Initial Solution Procedure
In the Nicolson-Ross-Weir procedure (5, 7) the equations for the scattering parameters are combined to allow the system of equations to decouple. This decoupling yields an explicit equation for the permittivity and permeability as a function of the scattering parameters on a point-by-point basis. This solution procedure is the basis of the commonly used techniques for obtaining permittivity and permeability. Unfortunately, these equations are not well-behaved for low-loss materials at frequencies that correspond to integer multiples of one half wavelength in the sample.
The Nicolson-Ross-Weir procedure that is implemented provides a good initial approximation to AuAz, B\ and Bz. All other parameters are initialized to zero. The estimated values for permittivity and permeability are determined on a point-bypoint basis by frequency. The corresponding scattering parameters are computed with these values and then compared to the observed values. The computed values for \jut and e.t that provide the closest agreement between the observed and predicted 5-parameter data are used as the initial values to the regression model. ' The justification of new terms added to each of Eqs. 2-5 is based solely on empirical evidence found when solving for the primary mode.

ODRPACK Orthogonal Distance Regression Package
Briefly, ODRPACK [2] is an implementation of a trust region Levenberg-Marquardt algorithm. This type of trust region approach adaptively determines the region in which the linear approximation closely resembles the nonlinear model. The procedure allows both an ordinary least squares model, in which the uncertainties are assumed to be only in the dependent variable, and, an orthogonal distance regression model, where uncertainties also exist in the independant variables.
First order derivatives for the Jacobian matrices can be numerically approximated (finite difference approximation), or can be user-supplied analytical derivatives. The procedure performs automatic scaling of the variables if necessary, as well as determination the precision of the model in terms of machine precision. The ODRPACK includes many other features that assist the user in the modeling process. The model automatically determines the number of digits in the model, checks analytical derivatives provided by the user, and automatically selects many of the input parameters for the user, if desired.
Iterations are stopped in ODRPACK when any one of three stopping criteria is met. Two of these indicate that the iterations have converged to a solution. Sum-of-squares convergence indicates that the change in sum-of-squares observational uncertainty is sufficiently small. Parameter convergence indicates that the change in the estimated parameters is sufficiently small. The third stopping criterion is a limit on the number of iterations.

Initial Conditions
Many of the input options for ORDPACK can be set to their default values, as is done in this model. The most significant input parameters for modeling permittivity and permeability are the initial values for A\, A3, Bi and B3. Sensitivity to the initial solution for these parameters is discussed below and the selection of initial settings is covered in Sec. 3.4. When higher modes are included, the solution for the primary mode is used as the initial guess for the higher mode model. For all of the parameters that define the physical model except those for ^lR and ER mentioned above, the standard laboratory values are used. All additional parameters are initialized to zero.

Sample Characteristics
A total of seven samples were modeled to determine permittivity and permeability with ANA twoport 5-parameter data. The sample characteristics appear in Table 1. Figure 2 contains selected plots of the observed ANA 5-parameter data versus frequency for sample 6.

Numerical Results
The first results reported are those for the primary mode model for all seven samples described In Sec. 4.3. The first set of plots (Fig. 3) include the observed 5-parameter data (dots) for the YIG sample (sample 4) from the ANA overlaid by the predicted data (line) found by the model. The corresponding residual plots^ for this sample appear in Fig. 4.
For the first four samples the predicted and observed data are nearly identical. The residual plot for cross-linked polystyrene (sample 1) as shown in Fig. 5 reveals the systematic uncertainty due to the ANA. Fig. 6 illustrates a 5-parameter primary mode solution for barium titanate mix 1 with a high dielectric constant and its corresponding residual plot. For three of the samples with high dielectric constants, the new model produces responses for the primary and one of the higher modes. See Fig.  7 for the 5-parameter real and imaginary components of 521 for the sample exhibited in Fig. 6 with the model for higher order modes.
For samples 5, 6, and 7 the predicted 5-parameter data provide realistic primary mode responses although the residual plots for these samples indicate that higher modes are present. With the higher order model described in Sec. 3.3 and the new solution found, the problem is resolved for the last three samples. The new solution specifies the additional parameters /3i,..., )84. Fig. 8 shows the results from the higher order model for sample 5. Additional work in this area suggests that this model is only an approximation.

ANA Systematic Uncertainty
The difference of the predicted 5-parameter and the observed values for samples with low dielectric values revealed systematic uncertainty. (For the high dielectric samples, other sources of uncertainty including higher mode responses dominate the ANA-induced uncertainty.) Additional tests revealed that the source of the uncertainty is not related to the material tested in the waveguide. In fact, uncertainties in the 5-parameter data for the cross-linked polystyrene sample closely resemble the 5-parameter data for an empty waveguide. Attempts to further identify and remove the systematic uncertainty are underway. Figure 5 contains the 5n data for an empty waveguide and also the residual plot for the cross-linked polystyrene sample.

Permittivity and Permeability Estimates
The estimates of permeability and permittivity for the various samples are defined by the parame-ters^!, .. .A4,Bi,... BA as a function of frequency. Several of the samples are magnetic (see Table 1). For these samples the permeability was allowed to vaty, with ODRPACK required to determine the (1,0) value. Slight variations in the value of permeability are apparent although the deviation from the true value is small. The solution procedure also provides the standard deviation and confidence intervals for each of the estimated parameters.

Robustness of the ODRPACK Procedure
The robustness of a mathematical procedure is its ability to find a locally optimal solution from a varied initial condition. The robustness of the entire permeability and permittivity procedure depends on the robustness of the ODRPACK procedure and, more significantly, the robustness of the math-ematical model. The existence of alternative optima in the mathematical model can limit the range of the initial conditions to produce a particular solution. In addition, singularities in the absence of alternative optima may force the solution procedure to fail to determine directions of improvement and cause a premature termination.
For the samples in this study, variability in the robustness of the procedure depended on the sample studied. For the materials with small dielectric constant, the procedure readily determined the correct solution for a variety of initial conditions. For materials with higher dielectric constant, the procedure often converged quickly, although the existence of alternative optima in the mathematical model often required the repeated solution with varied initial conditions before an acceptable solution was found.
In particular, for the cross-linked polystyrene sample, initial values for fiR and CR were set to (1,0) and (1,0) and were constant over the entire frequency range. The solution (1,0), (2.517,0.0018) was found after 50 iterations of the solution procedure. Initial conditions of (1,0), and (4.0,0.01) failed to determine a solution, while altering the imaginary part of €R to IxlO""* resulted in the desired final solution once again. For material with high dielectric constants, the range was smaller. For example for the barium titanate mk 1, initial conditions of (1.0,0.0) and (200.0,1.0) produced the converged value of (1,0) and (269.0,1.70) in 59 iterations.
One should note that the solution procedure can change the value of the estimated parameters by large amounts in the early iterations. Hence it is possible that for cross-linked polystyrene an initial solution farther from the desired solution may in fact find a solution merely because of the solution path taken by the procedure.

Conclusions and Future Directions
The nonlinear optimization procedure using twoport scattering parameters determined permittivity and permeability for a large number of samples (see Fig. 9). The added capability that permits variations in certain input parameters provides a mechanism to adjust for measurement uncertainties. The extension of the model to include higher order mode responses is quite useful for high dielectric constant materials.