Novel Application of Fast Simulated Annealing Method in Brushless Motor Drive (BLMD) Dynamical Parameter Identification for Electric Vehicle Propulsion

Abstract

Permanent magnet brushless motor drives (BLMD) are extensively used in electric vehicle (EV) propulsion systems because of their high power and torque to weight ratio, virtually maintenance free operation with precision control of torque, speed and position. An accurate dynamical parameter identification strategy is an essential feature in the adaptive control of such BLMD-EV systems where sensorless current feedback is employed for reliable torque control, with multi-modal penalty cost surfaces, in EV high performance tracking and target ranging. Application of the classical Powell Conjugate Direction optimization method is first discussed and its inaccuracy in dynamical parameter identification is illustrated for multimodal cost surfaces. This is used for comparison with the more accurate Fast Simulated Annealing/Diffusion (FSD) method, presented here, in terms of the returned parameter estimates. Details of the FSD development and application to the BLMD parameter estimation problem based on the minimum quantized parameter step sizes from noise considerations are provided. The accuracy of global parameter convergence estimates returned, cost function evaluation and the algorithm run time are presented. Validation of the FSD identification strategy is provided by excellent correlation of BLMD model simulation trace coherence with experimental test data at the optimal estimates and from cost surface simulation.

1. Introduction

High performance permanent magnet brushless motor drive (BLMD) systems are now widely used in electric vehicle (EV) propulsion [1, 2, 3], because of their higher power factor and efficiency, and are central to modern industrial automation [4] in such scenarios as aerospace systems control and maneuverability, numerical control (NC) machine tools and robotics. The benefits accruing [1, 3, 5, 6] from the application of such servodrives over other electric motor systems are higher power and better torque to weight ratio, a considerable saving of energy and higher precision control of torque, speed and position which promote better electric vehicle propulsion performance and optimal EV target ranging along with automation and control. This is due largely to the high torque-to-weight ratio and compactness of permanent magnet (PM) drives, lower heat dissipation and the virtually maintenance free operation of brushless motors in inaccessible locations when compared to conventional DC & AC electric motors. The controllers of these machine drives, which incorporate wide bandwidth speed and torque control loops, are adaptively tuned to meet the essential requirements of system robustness, high tracking performance and EV range acquisition without overstressing the hardware components [3, 7, 8]. An essential feature of the adaptive optimal control procedure in such BLMD-EV systems is the accurate dynamical parameter identification strategy used in conjunction with an accurate BLMD model [9, 10, 11, 12], where sensorless feedback current control loops are employed for reliable torque control during EV propulsion, with multi-modal penalty cost surfaces.

This book chapter, which is divided into three separate but interrelated sections, concerns BLMD dynamical parameter identification in which the choice of target data used in the Minimum Squared Error (MSE) objective function formulation has a significant effect on cost surface topology and selectivity in the vicinity of the global minimum [13, 14]. The penalty cost function selection in the identification of the motor drive dynamics impacts directly on the type of parameter search strategy to be adopted where embedded cost surface multiminima are concerned, which is discussed below, in terms of the accuracy of the returned parameter estimates. A short review of classical identification techniques is provided initially with an explanation of global convergence failure due to local minimum trapping over a multiminima cost surface based on BLMD step response feedback current (FC) target data. This difficulty with classical optimization methods highlights the need for an effective search strategy of parameter space with an in-built adaptive jump mechanism which can facilitate escape from possible local minimum capture and guarantee eventual global convergence. Such identification search features are provided by statistical optimization methods based on a Simulated Annealing (SA) kernel [15, 16, 17, 18, 19, 20, 21] which has an adjustable pseudo-temperature parameter that controls the jump related magnitude of the inherent random fluctuations. The Fast Simulated Diffusion algorithm [22, 23, 24], which is a more efficient version of SA that exploits the use of the classical gradient search technique at low pseudo-temperatures, surmounts the above obstacles and can be deployed as an effective global parameter extraction tool for BLMD identification.

In the first section the classical application of the Powell Conjugate Direction (PCD) optimization technique of motor parameter extraction over a two dimensional shaft velocity cost surface, which has a parabolic wedge shaped topography with a ‘line minmum’ stationary region [13, 14], is discussed. This will be used later for comparison purposes with the more successful Fast Simulated Diffusion (FSD) method, which is presented in the second section, to show that the FSD method is better in terms of the accuracy of returned FSD parameter estimates. It will also be shown that the PCD method does not converge to the global minimum and that the returned parameter estimates are less than satisfactory in the approximation of the optimal parameter vector when contrasted with the FSD method deployed over a multiminima FC cost surface. The convergence inaccuracy of the PCD method lays the groundwork for the introduction of a more accurate and efficient parameter identification search strategy based on the FSD method over a multiminima cost surface that is more selective at the global minimum.

In the second section the application of the FSD optimization search technique as an extension of the SA method is presented, with a multiminima objective function based on step response FC target data, for motor dynamical parameter extraction. Details of the FSD development and integration into the BLMD parameter estimation problem based on the minimum quantized step sizes from noise considerations are provided. The results of global parameter convergence estimates returned, which are very accurate, including the number of iterations, cost function evaluation and the algorithm run time are presented. Validation of the FSD identification strategy is provided by excellent correlation of BLMD model simulation trace coherence, at the optimal parameter estimates, with experimental test data and from cost surface simulation.

In the third section the development details of a novel modified form of the FSD algorithm are presented for fast accurate parameter identification [22, 23, 25, 26, 27, 28]. This will be based on the heuristics of the FC sinc-like surface topography in an effort to reduce the computational cost and search time to acquire global optimality. The beneficial effect of the incorporation of an approximated gradient search, along with the gathered cost statistics, at the end of each anneal step in the cooling schedule during the high temperature random search phase is discussed. Furthermore the condition, pertaining to the occurrence of random search trapping within the capture cross section of the cost surface containing the global extremum, for elimination of the reheat and thermal condensation phases of the FSD method is explained. Details of a set of tests to establish modified FSD (MFSD) performance in terms of convergence accuracy accompanied with a substantial reduction in search time, for three known cases of motor shaft load inertia (J) in the parameter identification process, are furnished with remote initialization in the vicinity of a potential local minimum trap. An error analysis of the returned parameter estimates is provided, which are almost identical to actual J values, and a comparison is made with alternative global estimates obtained from cost surface simulation to validate the modified FSD search technique. A theoretical analysis of the effect of FC data training record length on the probability of global minimum capture and capture cross sectional ‘area’ is provided along with a discussion of the impact of data record length on FC cost surface selectivity for SI purposes.

2. Motivation

The benefits accruing from the use of an embedded BLMD system in high performance adjustable speed drive applications has resulted in a need for an accurate physical model [12, 29] for the purpose of parameter identification of the drive dynamics. Allied to this need is a general requirement for an accurate and efficient search strategy of parameter space in the design of optimal drive controllers [30, 31, 32] where system identification is an implicit feature during online operation. This is necessary for PID auto-tuning of wide bandwidth current loops in torque control mode to speed up embedded BLMD commissioning and facilitate control optimization through regular retuning. The identification process is dependent on an accurate model of the nonlinear electromechanical system [33] which includes the pulse width modulated (PWM) inverter with power transistor turnon delay to avoid current shoot-through.

The extraction of the drive dynamical parameters generally relies on the minimization of some quantitative measure of error cost in terms of the goodness-of-fit, based on the MSE norm [13], between the observed motor drive output experimental test data and its model equivalent. The presence of multiminima in the MSE penalty function, however, results in a large spread of parameter estimates about the global minimum with model accuracy and subsequent controller design performance very dependent on the minimization technique adopted and the initial search point chosen. The existence of a noisy cost function, resulting in ‘false’ local minima proliferation in the stationary region containing the global extremum [13, 28], depends on the numerical accuracy with which the PWM delayed inverter switching instants are resolved in the model simulation [13]. Furthermore the plurality of genuine local minima is governed by the choice of data training record used in the objective function formulation which in the case of step response feedback current (FC) has a sinc-like topography [13]. The use of a step input (i/p) as a test stimulus is motivated by the fact that in normal industrial applications in the online mode input command changes are generally sudden and step-like and are sufficient for persistent excitation of the BLMD system. The accompanying transients in the observed variables can then be used effectively for parameter identification of motor shaft viscous damping factor and inertia changes during normal operation.

Experimental data training sets, of the observed variables [13] for various known shaft load inertia, can be used during the parameter identification procedure with the corresponding BLMD model simulation runs to establish a parameter mean squares error cost surface as the objective function to be minimized. The deployment of step response winding current feedback as a target function is found to be particularly beneficial as it exhibits the frequency modulated (FM) characteristic of a constant amplitude swept frequency sinusoid. The FC response also has better overall cost surface selectivity, which improves with observed data record length, by comparison with shaft velocity target information [13]. The corresponding simulated step response, upon parameter convergence, will be shown to have both frequency and phase coherence which attests to the accuracy of the identification methodology. Furthermore the selection of an FC multiminima cost function as a suitable choice for investigation in motor parameter extraction is motivated by the scenario, as a secondary consideration, that winding current flow information will be the only feedback signal available for sensorless motor control [34] applications.

The step response FC cost function exhibits an apparently smooth continuous one dimensional surface when plotted against either inertia or the damping factor as the free parameter to be identified [13] for coarse step size simulation. The local minimisers appear to be symmetrically disposed about the global minimum in accordance with a sinc function profile. The attendant unimodal velocity cost function appears to be parabolic in the inertia and friction parameter variables [13]. When high resolution of the motor shaft inertia extimates is required, in simulation trials, during system identification both cost functions appear to possess a granulated response surface thus rendering uncertainty in the parameter extraction process [13]. However this difficulty in parameter convergence is somewhat tempered by the observation that the noise floor in the error response surface is eclipsed by the residual error magnitude at the global minimum, which ameliorates the parameter search process, and is alleviated by the adoption of parameter quantization as discussed in [13, 28]. The response surface ‘noisiness’ arises primarily from the delay nonlinearity in the PWM current loops of the BLMD model structure and the accuracy with which the PWM crossover times are determined with the subsequent timing of the delayed inverter firing signals [9, 10]. The presence of ‘genuine’ local minimisers in the FC response cost surface is manifested through interference by the relative phasing of the swept frequency motor current sinusoids with target data at different inertia values [13, 28].

3. Choice of parameter identification methods

Two methods of parameter identification, which are based on the PCD and FSD optimization techniques and linked with the shape of the error response surface, will be investigated in sections 4 and 7 in the extraction process of known motor shaft inertia from data training records. The first, rooted in classical optimization techniques [35] where derivative information is not required [36, 37], relies on the application of Powell’s conjugate direction search to the concave velocity response surface with single parameter variation [13]. This classical parameter extraction technique is superior to and much more efficient [35], due to the orthogonality of its conjugate directions of search, than other direct search methods [38, 39, 40] such as the Simplex method [41] and the method of Hooke and Jeeves [42, 43] which are expensive in CPU time and slow to converge. Other more efficient classical optimization techniques [35, 44, 45, 46], such as the Polak-Riebere conjugate gradient method or the Newton-like BFGS method [47, 48] or the hybrid Levenberg–Marquardt method [49, 50, 51] are not considered in this case because of cost surface noisiness and the need to calculate partial derivatives which is difficult with the presence of ‘noise’ related PWM computational uncertainty. These alternative methods are computationally expensive in BLMD simulation time, where the evaluation of the gradient vector and Hessian matrix [35] are concerned, and are difficult to apply in any case because of the model complexity of the complete motor drive system with PWM inverter delay operation without resorting to difference equation approximation of the derivatives [35, 52]. Furthermore partial derivative evaluation is suspect in the presence of computational ‘noise’, inherent in the cost function, for an infinitesimal change in the relevant parameters. This can result in an erratic hill-descent, associated with the conjugate gradient direction search, over the response surface and entrapment of the minimization procedure in a false minimum in the noise grained incline of the bowl shaped velocity cost function [13, 28]. All classical minimization procedures are well known to have difficulty with the FC cost surface topography illustrated in [13, 28] because they are easily trapped in one of the embedded local minima and thus fail to converge to the optimal parameter set. This problem is accentuated by initialization of the search process remote from the global minimum, due to uncertainty, in a region where there are local minima. Classical hill-decent methods that rely on following the cost gradient are easily captured at local minima in this instance.

The second identification method is based on the statistical search technique of simulated diffusion (SD) [24, 53], as an adjunct of simulated thermal annealing (SA) [15, 16, 17] used in combinatorial optimization [25], and has to be deployed to acquire the global minimizer of the FC multiminima cost surface. The SD technique in contrast to the PCD method has been successfully applied to the multioptimal problem of MOSFET parameter extraction [19, 24] and simulated annealing in circuit placement tasks [20, 21], with a cost function having a fractal landscape, in the VLSI domain. Both the PCD and FSD methods of parameter extraction presented compare favorably in terms of the returned shaft load inertia, for short data records with initialization close to the global minimum, which enhances confidence in the model of the BLMD servodrive system as well as in the performance of both identification strategies. However the sensitivity of the FC objective function to inertia parameter changes is an order of magnitude greater than that of its shaft velocity counterpart [13, 28]. This results in better selectivity and more accurate convergence in the case of the SD method whereas Powell’s method is more susceptible to capture in a false minimum due to cost surface noise in the vicinity of the global minimum. Furthermore with longer data records the selectivity of the FC cost function increases to a maximum, when the motor has reached full speed, accompanied with ‘genuine’ local minima proliferation. The reverse effect is manifested as a flattening of the response surface in the neighborhood of the global minimum, with the resulting minimization procedure susceptible to remote trapping in noisy local minima, for the shaft velocity cost function. In the former procedure relative phase information can be used effectively for estimating parameter variation at motor speed saturation while the benefits of the initial speed response transient are attenuated in the latter case. The FSD search technique is verified, with initialization far from the global minimum, for known motor shaft inertia. The accuracy of the returned FSD estimates is also checked against the cost surface simulation. Convergence details and comparisons of both identification methods in terms of iteration count, objective function evaluations or motor simulation runs and CPU time are presented.

4. Application of Powell’s identification method to velocity cost surface

The choice of Powell’s Conjugate Direction Set method [35, 49, 54, 55] of unconstrained optimization as a method of BLMD dynamical parameter identification within the indicated parameter tolerance bounds was motivated by topographical considerations of the shaft velocity objective function E_ω(X) and ingrained cost surface noisiness. The penalty cost function concerned appears to have a ‘line minimum stationary region’ predominantly in the B parameter direction as elucidated in [13]. Consequently convergence difficulties can arise from application of any of the steepest decent conjugate gradient methods of parameter identification, such as the Polack-Ribiere technique [35], with this type of cost function. This problem results from the requirement that

$$\text{∇}{E}_{\mathrm{\omega }}{\left({\mathit{X}}^{\left(k+1\right)}\right)}^{\mathrm{T}}{\mathit{s}}^{\left(k\right)}<0$$E1

for hill decent and approaches zero as the global minimum is reached for the possible gradient search directions ∇E_ω(X_j^(k+1)) indicated in Figure 1.

Figure 1.

Global minima multiplicity.

The component of the gradient ∇E_ω(X^(k+1)), which is tangential to the cost contour in this instance, vanishes as shown along the k^th iterate search direction s^(k) at the points X_j^(k+1) for j = 1,2, etc. This condition results in a multiplicity of search directions and possible global minimum values X^(k+2) along the ‘line minimum’ of the wedge shaped syncline of the cost velocity surface as illustrated in Figures 2–4. This difficulty is partially borne out in [13] where the fitted response quadratic model parameters were re-evaluated at successive iteration points to improve the global convergence estimate. Furthermore the presence of point-like singularities due to cost surface noisiness, with ‘false’ local minima proliferation in the neighborhood of the global minimizer [13], results in a discontinuous cost function with consequent difficulties with derivative calculations. Also the derivative computation burden increases with expensive BLMD cost function evaluation where lengthy simulation times are concerned with a small time step Δt for accurate resolution of PWM edge transitions. The steepest decent gradient search technique can also suffer from oscillatory behavior and poor convergence results due to numerical round-off effects as reported by Fletcher [35].

Figure 2.

NSL velocity cost surface.

Figure 3.

MSL velocity cost surface.

Figure 4.

LSL velocity cost surface.

The PCD method, which relies on cost function evaluations only without derivative information, is cyclically deployed instead as a line search algorithm in parameter space beginning at X⁽¹⁾. This algorithm progresses along N mutually orthogonal directions conforming to the dimensionally of parameter space until sufficient accuracy of the global minimum estimate has been attained according to some convergence stopping criterion. With this method the iterate moves gradually towards the neighborhood of the global optimum X_opt via inexact line searches initially and then rapidly converges to the stationary point itself. The iteration process is terminated upon some user supplied convergence test of the form [49].

$$2\left\{(E_{\mathrm{\omega }}\left({\mathit{X}}^{\left(1\right)}\right)-E_{\mathrm{\omega }}\left({\mathit{X}}^{\left(N+1\right)}\right)\right\}/\left\{E_{\mathrm{\omega }}\left({\mathit{X}}^{\left(N+1\right)}\right)+E_{\mathrm{\omega }}\left({\mathit{X}}^{\left(1\right)}\right)\right\}<\mathrm{\epsilon }$$E2

being satisfied with a gradual reduction in the cost E(X) towards an accumulation point ${\widehat{\mathbf{X}}}_{\mathit{opt}}$ which approximates the global minimum X_opt within the specified error bound ε. Furthermore the application of a specific threshold step size δX^L [13, 28] for parameter space quantization ameliorates the difficulty with response surface noisiness [13, 28]. This methodology results in inexpensive line searching during parameter extraction with a reasonable degree of convergence accuracy maintained. In the application of the PCD algorithm a quadratic model is used to approximate the two dimensional MSE objective function, pertaining to the shaft velocity in terms of the motor dynamical parameters [13], so that a prediction of the location of the local minimum can be made. This localized response surface modeling technique guarantees second order convergence [35] and is a very suitable candidate for the velocity cost function given the parabolic nature of its topography with respect to the inertia parameter which is the most likely to vary in high performance drive applications. The applied PCD method is based on the property of quadratic termination of the approximation model, with Hessian $\widehat{\mathbf{G}}$>0, at the global minimizer. This model strategy, which is similar to the normal form in [13, 28], admits to the existence of at most N line searches {s⁽¹⁾, s⁽²⁾, …, s^(N)} in N-dimensional parameter space X for global convergence along independent mutually conjugate directions [35] such that

$${\mathbf{s}}^{{\left(i\right)}^T}\widehat{\mathbf{G}}{\mathbf{s}}^{\left(j\right)}=\left\{\begin{array}{c}1\forall i=j\\ 0\forall i\ne j\end{array}\right.$$E3

5. Application of PCD method to BLMD parameter extraction

The PCD method, which begins with a user supplied estimate X⁽¹⁾ and a set S of directions based on the column vectors of any orthogonal matrix such as the identity I with

$$\mathit{S}=\left\{{\mathbf{s}}^{\left(\mathrm{k}\right)}|{\mathbf{s}}^{\left(\mathrm{k}\right)}={\mathit{e}}^{\left(\mathrm{k}\right)}={\left[0,0,\dots {\mathrm{,0,1}}^{\left(\mathrm{k}\right)},0,\dots ,0\right]}^{\mathrm{T}},\forall k=1,\dots ,N\right\}$$E4

initially, proceeds according to a cyclic N dimensional search of parameter space as per the flowchart in Figure 5 until some stopping criterion based on sufficient cost reduction is satisfied. The basic structure of the algorithm can be summarized as a line minimum search routine due to Brent [49, 56] for the k^th iterative search as

• Determine the search direction s^(k) in discretized parameter space which is conjugate to all previous excursions into the quantized parameter lattice.

• Evaluate α^(k) so as to minimize E_ω(X^(k)+α^(k)s^(k)) with respect to α in discrete space X based on a quadratic model approximation of E_ω(X) in this direction and proceed to set

Figure 5.

Flowchart of Powell’s conjugate direction set optimization technique.

$${\mathit{X}}^{\left(k+1\right)}={\mathit{X}}^{\left(k\right)}+{\mathrm{\alpha }}^{\left(k\right)}{\mathit{s}}^{\left(k\right)}$$E5

with resolved accuracy given by the worst case quantized step size ±δX^L in [13, 28]. In the above line minimum iterative search procedure along a particular direction s^(k) the minimizer of the cost function is crudely bracketed initially by quadratic polynomial extrapolation. The cost function minimum estimate X^(k+1) is then obtained by successive approximation within the bracketed interval via interpolation from the fitted quadratic model Q(X) to the triplet of points {X₁, X₂, X₃}, where the response surface is concave, as

$${\mathbf{X}}^{\left(k+1\right)}={\mathbf{X}}_m={\mathbf{X}}_2-\frac{1}{2}\cdot \frac{{\left({\mathbf{X}}_2-{\mathbf{X}}_1\right)}^2\left[E_{\omega }\left({\mathbf{X}}_2\right)-E_{\omega }\left({\mathbf{X}}_3\right)\right]-{\left({\mathbf{X}}_2-{\mathbf{X}}_3\right)}^2\left[E_{\omega }\left({\mathbf{X}}_2\right)-E_{\omega }\left({\mathbf{X}}_1\right)\right]}{\left({\mathbf{X}}_2-{\mathbf{X}}_1\right)\left[E_{\omega }\left({\mathbf{X}}_2\right)-E_{\omega }\left({\mathbf{X}}_3\right)\right]-\left({\mathbf{X}}_2-{\mathbf{X}}_3\right)\left[E_{\omega }\left({\mathbf{X}}_2\right)-E_{\omega }\left({\mathbf{X}}_1\right)\right]}.$$E6

This minimum estimation procedure is assisted, if iterative progress towards quadratic termination is stalled in a non convergent limit cycle, by conducting a golden section search [56] of the bounded interval through contraction in terms of the golden mean g_m (∼0.382) as

$${\mathbf{X}}_m=\left\{\begin{array}{c}\left({\mathbf{X}}_3-{\mathbf{X}}_2\right)g_m;\text{If}{\mathbf{X}}_2<\left({\mathbf{X}}_1+{\mathbf{X}}_3\right)/2\\ \left({\mathbf{X}}_2-{\mathbf{X}}_1\right)g_m;\text{If}{\mathbf{X}}_2>\left({\mathbf{X}}_1+{\mathbf{X}}_3\right)/2\end{array}\right.$$E7

to trap the least MSE cost E_ω(X^(k+1)). During a single cycle iteration of the PCD algorithm the single largest cost decrease ΔE_ω(X) along a particular search direction s^(dm) is monitored. This cost reduction metric is then used to determine whether or not the set S of conjugate directions needs to updated before commencing the next Powell (J+1)^th iteration cycle. The decision to include a new search direction, by replacing that in S along which the largest cost decrease ΔE_ω was observed, is made by first evaluating the cost at the extension point

$${\mathit{X}}^{\mathit{EXT}}={\mathit{X}}^{\left(k+1\right)}+{\mathit{s}}^{\mathit{AV}}$$E8

along the proposed average direction

$${\mathit{s}}^{\mathit{AV}}={\mathit{X}}^{\left(k+1\right)}-{\mathit{X}}^{\left(k\right)}$$E9

to avoid linear dependence windup and resultant loss of conjugacy. This cost is then compared with that at the initial search point X^(S) along with the differential comparison

$$\mathrm{\Delta }E=\left[(E^{\left(S\right)}-E_{\omega }\left(X^{\left(k+1\right)}\right)\right]-\mathrm{\Delta }{E}_{\omega }$$E10

to determine if any reduction is achievable along s^AV for the J^th iteration, which is regarded as the average path traversed over all possible N directions in S, other than that via s^(dm). Furthermore the curvature estimate at X^(k+1) in the direction s^AV, given by

$$\kappa =\left\{E_{\omega }\left(X^{\left(S\right)}\right)-2E_{\omega }\left(X^{\left(k+1\right)}\right)+E_{\omega }\left(X^{\mathit{EXT}}\right)\right\}$$E11

is also checked for possible minimum convergence of the iterative search along with Eq. (10) using the flowchart test condition in Figure 5 as

$$2\kappa \cdot \mathrm{\Delta }{E}^2\ge \mathrm{\Delta }{E}_{\omega }{\left\{E_{\omega }\left(X^{\left(S\right)}\right)-E_{\omega }\left(X^{\mathit{EXT}}\right)\right\}}^2.$$E12

If the test condition Eq. (12) is true then the old set of directions is retained for the next Powell iteration either because cost reduction is already exhausted or no significant reduction is observed in any particular direction or a substantial second derivative exists along s^AV indicating minimum convergence of the search estimate. Conversely if Eq. (12) is false then the direction set S is updated by replacing s^(dm) with s^AV before embarking on the next pass of the Powell algorithm which is usually restricted to some maximum user supplied iteration count N_max of typically 200.

The PCD method was applied to the identification of three known values of BLMD shaft load inertia with experimental shaft velocity target data deployed in the MSE objective function E_ω(X) formulation displayed in Figures 2–4. An insight into the progress of this optimization method towards global optimality in the extraction of the dynamical J and B parameters with initialization at the tolerance band edge can be obtained, for example with zero shaft inertial load conditions, from Figure 2 over the generated velocity cost surface. After one cycle of the Powell method the iterate has reached the ‘line minimum’ of the cost surface syncline by completing alternate line searches along the J and B parameter directions. Quadratic convergence of the PCD method thereafter is relatively ‘slow’ in that a further three iteration steps of the PCD algorithm are required to reach a limit point estimate ${\widehat{\mathbf{X}}}_{\mathit{opt}}$ of the global minimum. This is manifested as a zigzag search pattern over the stationary region of the response surface in the quantized J and B parameter directions along the ‘line minimum’, which is essentially in the B direction, where the cost function curvature is low with inadequate selectivity for global minimum reachability.

6. Parameter convergence results for the PCD method

The optimal estimates of the BLMD dynamics returned by the PCD method including convergence details for three known cases of shaft load inertia are summarized in Table 1. The penalty cost reduction sequence associated with the application of the Powell algorithm is displayed in Figures 6–8 for each of identified shaft inertial loads along with the cumulative number of bracketing intervals and cost function evaluations required at each iterative step. Substantial cost reduction, which can be attributed mostly to the adjustment of the inertial parameter, ceases after one iterative cycle of the PCD method in each case signifying the arrival of the parameter estimate in the stationary region along the ‘line minimum’ as in Figure 2.

Total Shaft Load Inertia; J (kg.cm^2)	Figure 2 - NSL; No Shaft Load; Jm = 3.0	Figure 3 - MSL; Medium Inertia; J = 12.303687	Figure 4 - LSL; Large Inertia; J = 20.822
Parameter Initialization: {Jinit = 0.82 J, Binit = 1.09Bm} with Bm = 2.14x10^−3 Nm/rad/sec
Initial Cost Einit	10.2481x10^−2	36.5942x10^−3	54.4345 x10^−3
Returned Parameter Estimates after first Iteration Cycle of PCD Method
$J_{\mathit{opt}}^{\left(1\right)}$ (kg.cm^2)	3.02659	12.0721	21.292
$B_{\mathit{opt}}^{\left(1\right)}$ (Nm.rad^−1.sec)	2.10405x10^−3	2.14214x10^−3	2.06596x10^−3
Optimal Parameter Estimates returned from Powell’s Conjugate Direction Method
${\widehat{J}}_{\mathit{opt}}$ (kg.cm^2)	3.168237	12.29243	21.51233
${\widehat{B}}_{\mathit{opt}}$ (Nm.rad^−1.sec)	1.913588x10^−3	1.989772x10^−3	1.95168x10^−3
Min. Cost ${\widehat{E}}_{\mathit{opt}}$	3.278108x10^−2	5.910952x10^−3	2.456056x10^−2
Total No. of PCD Iterations	4	4	2
No. of Cost Func.-Eval^s	99	121	56
Initial Bracketing	24	31	15
Line Searching	74	89	40
Av. No. of Func Eval^s/Iter^n	∼25	∼30	28
Total Time (sec)	3584	5234	2918
Average Iter^n Time (sec)	896	1308.5	1459
Time/MSE Cost Eval. (sec)	36.202	43.256	52.106

Table 1.

Returned BLMD dynamical parameter estimates via PCD algorithm.

Figure 6.

NSL iterative cost reduction.

Figure 7.

MSL iterative cost reduction.

Figure 8.

LSL iterative cost reduction.

This observation can be adduced from Table 1 where the listed parameter estimates after the first iteration, in each test case, are reasonably close to the eventual global minimizer estimates. Further gains in MSE reduction are expensive in this zone of the cost surface, due to poor selectivity as explained in [13], with the bulk of the PCD effort devoted towards improving global minimum convergence estimate of the returned parameter vectors. The computational runtime of Powell’s algorithm increases in accordance with the iterative count and the accumulated number of BLMD simulation related MSE evaluations for each of the inertial test cases displayed in Figures 9–11. The corresponding averaged BLMD simulation time increases almost in proportion with the inertial loading for a fixed experimental shaft velocity data capture V_ωr, normalized to about 10 equivalent machine cycles of motor FC transient as in Tables 2 and 3, with time decimation to 4095 sample points.

Figure 9.

NSL cumulative iteration time.

Figure 10.

MSL cumulative iteration time.

Figure 11.

LSL cumulative iteration time.

δJ^L = 15.7386x10^−3 kg.cm^2 ≡ 0.51%Jm; Jm = 3.086x10^−4 kg.m^2 (NSL - Rotor Inertia)	δB^L = 38.092x10^−6 Nm.rad^−1.sec ≡ 1.78%Bm; Bm = 2.14x10^−3 Nm.rad^−1.sec
Shaft Velocity Target Data Vωr; Data Sample Rate TS; Decimation Factor	No Shaft Load; 20 μs; 20	Medium Inertia; 40 μs; 40	Large Inertia; 49.6 μs; 50
No of Equivalent FC Cycles; for Computation Benchmarking	∼9.75	∼11.5	∼10.5
Simulation Grid Size in Eω; [ΔJ, ΔB] as per [13]	Figure 2; [2δJ^L, δB^L]	Figure 3; [4δJ^L, δB^L]	Figure 4; [7δJ^L, δB^L]
Surface Minimum Estimates; ${\overline{J}}_{\mathit{opt}}$ (kg.cm^2); ${\overline{B}}_{\mathit{opt}}$ (Nm.rad^−1.sec); ${\overline{E}}_{\mathit{opt}}$	3.184; 1.875 x 10^−3; 3.26 x 10^−2	12.38; 1.875 x 10^−3; 5.819 x 10^−3	22.36; 1.875 x 10^−3; 2.453 x 10^−2

Table 2.

Details of velocity cost surface generation with parameter quantization.

This capture restriction of the shaft velocity target data is due to the experimental constraints of the data acquisition system used, which limited the size of the data record acquired. These returned PCD statistics are also based on a fixed time step of 1 μs, in relation to BLMD model exercise in the MSE penalty cost formulation, for accurate realization of the PWM edge transitions and for computation benchmarking purposes in both the PCD and modified FSD methods of identification.

The shaft velocity data sets can also be deployed in the construction of MSE response surfaces E_ω over the two dimensional [J,B] quantized parameter manifold [13], with mesh size details given in Table 2, as a secondary means of parameter extraction for qualification of the PCD method and verification of the returned PCD optimal estimates in terms of accuracy. These cost surfaces, pertaining to the various motor shaft inertial loads, are illustrated in Figures 2–4 and appear to posses a parabolic ravine-like structure with an embedded elliptical stationary region stretched into a ‘line minimum’ shape [13]. These cost constructs can be used as the basis for an error analysis as given in Table 4 by which the accuracy of PCD estimates are gauged.

FC Target Data; No. of Cycles; Data Sample Rate Ts; Decimation Factor	No Shaft Load; ∼9.75; 20μS; 20	Medium Inertia; ∼11.5; 40μS; 40	Large Inertia; ∼10.5; 49.6μS; 50
Simulation Mesh Size; [ΔJ, ΔB]	$\left[2\delta J^L,\delta B^L\right]$	$\left[6\delta J^L,\delta B^L\right]$	$\left[10\delta J^L,\delta B^L\right]$
Surface Minimum; ${\overline{\mathrm{J}}}_{\mathrm{opt}}$ ($\mathrm{kg}.{\mathrm{cm}}^2$); ${\overline{B}}_{\mathit{opt}}$$\left(\mathrm{Nm}.{\mathrm{rad}}^{‐1}.sec\right)$	3.089; 1.921$\times {10}^{‐3}$	12.158; 1.921$\times {10}^{‐3}$	20.877; 1.921$\times {10}^{‐3}$
Min Cost ${\overline{E}}_{\mathit{opt}}$	4.922$\times {10}^{‐2}$	1.877$\times {10}^{‐2}$	1.236$\times {10}^{‐2}$

Table 3.

Quantized parameter FC response surface simulation details.

Motor Inertia Jopt; $\left(\frac{{\widehat{J}}_{\mathit{opt}}-\overline{J}}{\overline{J}}\right)\times 100\%$	No Shaft Load; −0.495%	Medium Shaft Inertia; −0.707%	Large Shaft Inertia; −3.79%
Damping Factor Bopt; $\left(\frac{{\widehat{B}}_{\mathit{opt}}-\overline{B}}{\overline{B}}\right)\times 100\%$	−2.058%	−6.121%	−4.09%
Global Cost Eopt; $\left(\frac{{\widehat{E}}_{\mathit{opt}}-\overline{E}}{\overline{E}}\right)\times 100\%$	0.555%	1.58%	0.125%

Table 4.

Error analysis for returned PCD parameter estimates (% relative error).

The percentage relative error in the returned PCD optimal estimate ${\widehat{J}}_{\mathit{opt}}$ appears to increase with shaft inertial load but with no apparent trend noticeable in the damping parameter estimate. The percentage error in the identified ${\widehat{B}}_{\mathit{opt}}$ parameter, however, exceeds that for the inertia in all three test cases which highlights the problem in accurately extracting the damping coefficient due to presence of a pronounced ‘line minimum’ of admissible friction values predominantly in the B parameter direction.

The opposite error pattern will be shown to occur in the extraction of the inertia estimate, coupled with lower relative error in the identification of the friction coefficient, upon application of the modified form of the FSD method with details given in Table 5. This contrast in error pattern for the returned FSD parameters is due to better discriminating features of the corrugated FC cost surface in terms of its selectivity, with greater curvature in the B parameter direction for increased data capture length, and lower threshold parameter step sizes as shown in Tables 7–9 in [13].

Motor Inertia Jopt; $\left(\frac{{\widehat{J}}_{\mathit{opt}}-\overline{J}}{\overline{J}}\right)\times 100\%$	No Shaft Load; −0.29%	Medium Shaft Inertia; −0.15%	Large Shaft Inertia; −0.19%
Damping Factor Bopt; $\left(\frac{{\widehat{B}}_{\mathit{opt}}-\overline{B}}{\overline{B}}\right)\times 100\%$	1.98%	4.16%	1.98%
Global Cost Eopt; $\left(\frac{{\widehat{E}}_{\mathit{opt}}-\overline{E}}{\overline{E}}\right)\times 100\%$	5.99%	17.37%	5.42%

Table 5.

Error analysis for returned FSD parameter estimates (% relative error).

Nominal Shaft Inertial; Load; Jnom ($\mathrm{kg}.{\mathrm{cm}}^2$)	No Shaft Load; (NSL - Rotor); Jm = 3.0	Medium Shaft Load; (MSL); 12.304	Large Shaft Load; (LSL); 20.822
Current Feedback; Ifa	Figure 12; 47.4%	Figure 13; 67.2%	Figure 14; 61.0%
Current command; Ida	Figure 15; 44.2%	Figure 16; 67.6%	Figure 17; 53.0%
Current Controller o/p; Vca	Figure 18; 18.1%	Figure 19; 38.8%	Figure 20; 20.4
Motor Shaft Velocity; Vωr	Figure 21; 93.7	Figure 22; 99.1%	Figure 23; 98.6%

Table 6.

BLMD simulation trace coherence correlation coefficient ρ.

Data Record Length Nd	4095	8000	12000	24000	32000
Surface Eω(${\overline{J}}_{\mathit{opt}}$) Curvature ${\kappa }_J^E$	3.122x10^7	2.843x10^7	2.036x10^7	1.026x10^7	7.692x10^3
Convergence Metric $N_J^{\delta }$	10.3	10.8	12.8	18	20.8
Best Parameter Resolution Possible with PCD Algorithm for δJ^L = 15.739x10^−3 kg.cm^2
Total Shaft Inertial Load J; ($\mathrm{kg}.{\mathrm{cm}}^2$)		Shaft Inertia; Jm = 3.0	Medium Inertia; 12.304	Large Inertia; 20.822
% parameter Resolution δJ^L/J		0.525%	0.128%	0.0756%

Table 7.

Cost surface selectivity measure – Figures 24,25 and 26.

Ifa Cost Surface Target Data [13]		Sample Size Ni = 200 ≈ 10%NC
Parameters x	Inertia Jm kg.cm^2	Damping Bm Nm/rad/sec
Nominal Values xm	3.0	2.14x10^−3
Parameter Bounds ±Δxm	±20%	±10%
Quantized Step Size δx^L	9.1125x10^−3	2.6712x10^−5
Population Size Nx	141	16
Results of Initial Exploratory Search of Parameter Space
Histograms	Figure 30	Figure 31
Class Interval Size	10δJ^L	δB^L
Minimum Value xmin	2.408	1.953
Maximum Value xmax	3.592	2.327
Actual Mean Value xmean	3.041	2.153
Theoretical Mean Value $\overline{x}$	3.0163	1.9934
% Error in Mean Value	0.82%	8%
Actual Standard Deviation ${\widehat{\sigma }}_x$	0.343	0.109
Theoretical Standard Deviation ${\sigma }_x$	0.518	0.750
Optimal Parameter Estimate ${\widehat{x}}_{\mathit{\text{iopt}}}$	3.06379	1.953
FSD Exploratory Search Results of FC - MSE Cost Surface
Initial Global Minimum Estimate ${\widehat{E}}_{\mathit{\text{iopt}}}$		5.473x10^−2
Maximum Value Emax		6.936x10^−1
Mean Value $\overline{E}$		3.6671x10^−1
Standard Deviation σ		1.6538x10^−1
Initial Temperature $T_i$ (= kσ for k = 10)		1.6538

Table 8.

Initial details of FSD exploratory search phase.

Data Record Length Nd	6280	9000	12000	24000	32000
Fitted Coefficient b0	3.708x10^9	7.053x10^9	9.495x10^9	1.93x10^10	2.398x10^10
Selectivity Measure $S_J^E$	3.379x10^3	6.427x10^3	8.652x10^3	1.759x10^4	2.186x10^4
Parameter Resolution Accuracy in FSD Estimation with δJ^L = 9.1125x10^−3 kg.cm^2
Total Shaft Inertial Load J (kg.cm^2)		Shaft Inertia; Jm = 3.0	Medium Inertia; 12.304	Large Inertia; 20.822
% parameter Resolution δJ^L/J		0.304%	0.074%	0.0438%

Table 9.

Cost surface selectivity measure.

The observed phase-a current demand, feedback and controller o/p waveforms, displayed in Figures 12–20 at critical internal nodes of the BLMD model [10], appear to be to be coherent with similar test data in the early phase of the transient response of an actual drive system to a unit torque demand step i/p. However there is an eventual loss of synchronism, with the evolution of the model step response towards steady state conditions, in the waveform comparison due to the impact of the estimated J/B dynamic time constant mismatch in relation to the intrinsic value τ_m of the actual motor drive system.

Figure 12.

NSL FC simulation.

Figure 13.

MSL FC simulation.

Figure 14.

LSL FC simulation.

Figure 15.

NSL current demand simulation.

Figure 16.

MSL current demand Simulⁿ.

Figure 17.

LSL current demand Simulⁿ.

Figure 18.

NSL - current controller o/p.

Figure 19.

MSL - current controller o/p.

Figure 20.

LSL - current controller o/p.

This is borne out by the low correlation measurement coefficients in Table 6, which gauges the degree of FM coherence, between the respective experimental and simulated trace responses. Further evidence of this mismatch, though small, can be visualized in the deviation of the shaft velocity characteristics depicted in Figures 21–23, despite the high correlation coefficient to the contrary, as the motor drive accelerates towards rated shaft speed. A possible explanation for the correlation discrepancy between waveform types may attributed to the application in this instance of shaft velocity target data in BLMD parameter extraction process which results in a good fit between shaft velocity waveforms based on returned estimates due to PCD cost minimization. The returned estimates, however, result in a less than satisfactory trace coherence measure of fit in the current related characteristics pertaining to the BLMD current loop operation in toque control mode. Furthermore the use of larger quantized parameter step sizes, in conjunction with shaft velocity target data in the PCD parameter extraction process, results in a greater spread of returned optimal estimates about the global extremum compared with the threshold values deployed with current feedback data in the FSD algorithm.

Figure 21.

ZSL - shaft velocity simulation.

Figure 22.

MSL - shaft velocity simulation.

Figure 23.

LSL - shaft velocity simulation.

The initial impact of the shaft velocity transient step response decays and is ultimately swamped with the onset of steady state conditions as maximum shaft speed is reached with lengthening data records. The quadratic shape of the response surface E_ω is de-emphasized with reduced curvature, as more velocity target data is accumulated for MSE cost formulation, resulting in a loss of selectivity at the global minimizer.

This flattening of the cost profile about X_opt is clearly evident in Figure 24 for response surface cross sections in the inertia parameter J, based on simulated BLMD model target data, with increased data record length. The resultant ill conditioned surface admits a multitude of possible global minimum estimates about X_opt. The variation in cost function selectivity can be obtained by fitting a quadratic polynomial approximation [13] to each of the response surface cross sections E_ω(J) in Figure 24 and estimating the curvature at the fitted parabolic vertex ${\overline{J}}_{\mathit{opt}}$ as

$${\kappa }_J^E={\left.d^{2}Q\left(j\right)/{\mathit{dJ}}^2\right|}_{J={\overline{J}}_{\mathit{opt}}}=2b_2$$E13

in Table 7 with

$$Q\left({\overline{J}}_{\mathit{opt}}\right)=0\text{and}{\left.\mathit{dQ}\left(J\right)/\mathit{dJ}\right|}_{J={\overline{J}}_{\mathit{opt}}}=0.$$E14

Figure 24.

Cost surface variation with data record length.

This curvature ${\kappa }_J^E$ variation, shown in Figure 25, exhibits a decreasing quasi-linear dependency with target data length and indicates as a consequence poor response surface selectivity and large convergence radius in the stationary region containing the global minimum. An estimate of the global convergence radius r_J of the J parameter stationary zone can be obtained in terms of the fixed cost noise estimate $\widehat{\sigma }$, due to inexact PWM simulation [13], and the curvature variation with different data training record sizes N_d via Eq. (13) as

$$r_J=\sqrt{\frac{8\widehat{\sigma }}{{\kappa }_J^E}}.$$E15

Figure 25.

Cost surface selectivity.

This radial width can be referenced to the threshold step size δJ^L in Table 2, in order to establish the scale of the bounded region of convergence, with metric variation for different data sequence lengths N_d given in Table 7 as

$$N_J^{\delta }=r_J/\delta J^L$$E16

and displayed in Figure 26. This characteristic illustrates clearly the emergence of a pattern of reduced cost selectivity, which is mirrored as a radial extension of the global minimum region, with increased target data capture. The opposite trend prevails in Figure 27 with the application of current feedback in the MSE penalty cost and is the main motivation for its use as a target function with a tighter bound, besides sensorless motor control issues, in the identification of the BLMD dynamics despite the FSD computational intensity. The use of a fixed threshold step size establishes the degree of accuracy possible in parameter resolution during BLMD system identification as shown in Tables 7 and 9. The relative percentage accuracy in the returned estimates improves with motor shaft inertial loading as tabulated with a greater resolution possible with the deployment of current feedback target data in MSE cost reduction. If a coarser step size δX^* is adopted, which may be tolerable at large inertial loads without adversely affecting the percentage resolution in the returned estimates, both the FSD and PCD methods of parameter extraction can proceed much faster with smaller computational burdens. This results from the reduced number of the feasible lattice points to be searched in the quantized parameter domain with small changes in the percentage accuracy. The required parameter accuracy can be user defined in such circumstances at the start of the identification search routine and encoded in the relevant step size δX^* as a measure of the desired coarseness of resolution.

Figure 26.

Global convergence variation.

Figure 27.

Motor shaft speed variation.

7. Description of FSD method of parameter identification

The simulated annealing technique used here as a global optimization algorithm for motor parameter identification [22, 23] is based on the Fast Simulated Diffusion (FSD) method proposed by Sakurai et al. [24]. This method, motivated by quantum mechanics (QM), is modeled on the diffusion $\nabla E\left(x_i\right)$ of a particle, with position trajectory x_i, across a barrier potential E(x_i) under the influence of a temperature controlled random driving force $\sqrt{2T}{\mathit{dw}}_i$ resulting in Brownian motion. The particle ensemble dynamics [53] can be described by the Ito stochastic differential equation with space time coordinates (x, t) as

$$\mathit{dx}=-\nabla E\left(x\right)\mathit{dt}+\sqrt{2T}\mathit{dw}.$$E17

The down-hill gradient term -∇E(x) in Eq. (17), which is effective at low temperature, has the tendency to minimize the particle potential $E\left(x\right)$ along the continuous “parameter” path x_i as

$$\frac{{\mathit{dx}}_i}{\mathit{dt}}=-\nabla E\left(x_i\right)$$E18

and thus improve the cost of the objective function E(x). The stochastic perturbation component of motion ${\mathit{dw}}_i$ along a particular trajectory x_i(t), which is an essential coherent energy interaction at high temperature T, imparts enough momentum to enable the particle to “hill climb” its way out of a potential well and thus avoid local minimum capture with

$${\mathit{dx}}_i=-\nabla E\left(x_i\right)\mathit{dt}+\sqrt{2T}{\mathit{dw}}_i.$$E19

If a proper cooling procedure [57] is implemented the global minimum potential E_min is reached after infinite time with a Gibb’s distribution having a limiting Boltzman probability density [58].

$$p\left(x\right)\propto exp\left(-\frac{E\left(x\right)}{T}\right)\mathrm{as}t\to \infty$$E20

which peaks with a value of one at the global minimizer as T approaches zero. The global minimizer of E(x) can be obtained according to Aluffi-Pentini et al. [53] from inspection of the asymptotic value of a numerically computed sample trajectory x(t) as time t → ∞ via numerical integration of Eq. (17) from x(0) at t = 0 as the stochastic temperature T is reduced very slowly to zero with time t.

The differential random process in Eq. (17) has been mapped into an algorithmic procedure for approximating the global minimizer [24] of multiminima objective functions with a large N-Dimensional parameter manifold. Notable applications of discrete SA type methodology include the time scheduling problem of a traveling salesman to N cities which is NP-complete [15] and to cell placement in integrated circuit design. The SA method, however, requires very large computing resources [49] for global minimum convergence in such problem solutions whereas the introduction of the gradient search component in Eq. (17) with an aggressive cooling schedule can results in significant gains in the reduction of CPU time. In general the true location in parameter space of the global minimum of a multiminima cost function cannot be guaranteed with certainty in practical applications of the SA algorithm but a good sub-optimal approximation can be obtained with reasonable computational effort which is an acceptable estimate ${\widehat{\mathbf{X}}}_{\mathit{opt}}$ of the optimal parameter vector X_opt. Similar comments pertain to the FSD global minimum approximation during BLMD parameter extraction over the FC cost surface which is described herein.

In BLMD model parameter optimization the particle energy E(x) is replaced by the FC cost function E_Ifa(X) and its positional x_i configuration space by the parameter set X bounded by the permissible tolerance band ± ΔX_m. The FSD algorithm relies on two modifications, instead of integrating Eq. (19) directly, to accelerate convergence and prune wayward or non profitable moves in parameter space. The first is based on the use of an accept/non-accept function rule after Metropolis et al. [59] where the probability of acceptance of the next move

$${\mathbf{X}}_{k+1}={\mathbf{X}}_k+\mathit{d}\mathbf{X}$$E21

in the parameter identification search is governed by the Boltzman’s distribution. If the next move yields a lower cost E(X_k+1) than the current value E(X_k) at ${\mathit{X}}_k$ then the new parameter configuration X_k+1 is accepted. Alternatively in an error increasing move a random number

$$R_n\in \left[0,1\right]$$E22

is chosen for arbitration in the move selection process and the transition probability

$$P_{k+1}=exp\left(\frac{-\left\{E\left({\mathit{X}}_{k+1}\right)-E\left({\mathit{X}}_k\right)\right\}}{T}\right)$$E23

is computed. If the resulting decision is

$$P_{k+1}>\mathit{RN}$$E24

then the displacement X_k+1 is accepted, as a controlled uphill step as part of the iterative improvement process for a better solution, otherwise it is rejected in which case X_k+1 has to be regenerated. The second modification involves the alternate application of the hill-descending gradient search term -∇E(X) with the temperature dependent stochastic term in Eq. (19) in the generation of the next move. This alternate application ensures limited hill-descent at high temperature than otherwise would be the case if the two terms were added together as in the conventional method where ineffective moves due to the stochastic term would negate gains made in the gradient search. Instead of calculating the direction of steepest descent -∇E(X), which is expensive for large dimensional parameter space, the alternative directional entity

E25

is used instead because its expected value approaches -∇E(X) in the long term where $\widehat{\mathbf{r}}$ is a unit vector along a randomly chosen axis in parameter X space. A quadratic fit is then employed if the cost function is concave along the randomly chosen axis and the minimum is estimated via the Newton–Raphson method. If the cost function happens to be convex in the chosen direction a small $\mathit{d}\mathit{X}$ is first used and then doubled up a prescribed number of times until $E\left({\mathit{X}}_k+\mathit{d}\mathit{X}\right)$ fails to decrease. This approach gives a crude but inexpensive estimate of the minimum.

8. Application of FSD method to motor parameter estimation

The FSD iterative process proceeds slowly in accordance with a pseudo temperature control parameter T_k, with the same units as the cost surface, from an initial high temperature value T_i. This annealing technique permits the generation of a population of parameter vector X adjustments at each anneal temperature step to simulate the effect of cooling tardiness for the parameter ensemble to reach a steady state configuration.

The cooling process also includes an equilibrium condition based on some minimum level of acceptance at each temperature step before progressing on to the next step. The temperature $T_k$ is reduced very slowly, according to a prescribed anneal curve [60] or exponentially [15], to avoid premature quenching and resultant parameter configuration trapping in a metastable state. The outer temperature loop decrementation sequence and parameter adjustment acceptances at each inner loop equilibrium temperature constitute the simulated annealing schedule. The temperature reduction should be rapid with an early detection of equilibrium for effective cooling without quenching. Condensation followed by termination of the anneal method occurs when the cost function value tends to reduce slowly and remain unchanged for several consecutive temperature steps according to some stopping criterion. At this stage the parameter configuration has frozen to an optimal arrangement in the neighborhood of the global minimum cost.

The flowchart of the complete FSD algorithm is presented in Figures 28 and 29 with details of the inner and outer temperature loops. The high temperature melting phase, used in initialization, and the cooling procedure for the FSD algorithm are adopted from Huang et al. [60]. The annealing start temperature T_i is statistically determined from the standard deviation σ of the cost function distribution over N_i sample points from an initial exploratory search of parameter space. During this hot phase the feasible parameter domain is uniformly sampled with the temperature T_i assumed high enough at ‘infinity’ $T_{\infty }$ such that all generated states are accepted. An adequate sample size is established from the tolerance bounds ±ΔX_m imposed, which define a hypercube of feasible lattice points in parameter space, with interstitial distance based on the quantized step sizes [13]. A sample size of 10% of the total number determined [13, 28] as

Figure 28.

Flowchart of fast simulated diffusion algorithm – Heating phase.

Figure 29.

Flowchart of fast simulated diffusion algorithm – Cooling phase.

$$N_C=N_B·N_J=\left(\frac{2\mathrm{\Delta }{B}_m}{\delta B^L}\right)\left(\frac{2\mathrm{\Delta }{J}_m}{\delta J^L}\right)=\left(16\right)·\left(141\right)=2256$$E26

respectively, gives a good search coverage of parameter space for tolerance bounds given in Table 8 with $N_i=200$. This is evident from the histograms of the near uniform parameter search distributions shown in Figures 30 and 31 with details in Table 8 for interval sizes of δB^L and 10δL^L. The uniformity of the exploratory search distribution can be checked from theoretical consideration of the population size N_x of quantized parameter values, for the given tolerance bounds, as

Figure 30.

Histogram of J search value.

Figure 31.

Histogram of B search value.

$$N_x=\left(\text{2Δ}{x}_m\right)/{\mathit{\delta x}}^L$$E27

The mean parameter estimate obtained from a random search of discretized parameter space, with sample space size N_x and uniform probability of occurrence

$$p_x=1/N_x,$$E28

is given by

$$\begin{array}{ll}\overline{x}& =\frac{1}{N_x}\sum _{j=1}^{N_x}{X}_i=\frac{1}{N_x}\sum _{j=1}^{N_x}\left[\left(x_m-\mathrm{\Delta }{x}_m\right)+\left(j-1\right)\delta x^L\right]\\ & =\left(\frac{N_x-1}{2N_x}\right)\left\{2\left(x_m-\mathrm{\Delta }{x}_m\right)+\left(N_x-1\right)\delta x^L\right\}\end{array}$$E29

This theoretical estimate compares favorably with those in Table 8 obtained from FSD simulation with low relative error percentages, which verifies the randomness quality of the initial search. The standard deviation ${\sigma }_x$ of the parameter search estimates can be likewise determined theoretically via Eqs. (28) and (29) from the variance ${\sigma }_x^2$ with

$$\begin{array}{l}{\sigma }_x^2=\overline{x^2}-{\overline{x}}^2=\frac{1}{N_x}\sum _{j=1}^{N_x}{x}_j^2-{\overline{x}}^2\\ ={\left(x_m-\mathrm{\Delta }{x}_m\right)}^2+\left(N_x-1\right)\cdot \left(x_m-\mathrm{\Delta }{x}_m\right)\delta x^L+\frac{{\left(\delta x^L\right)}^2}{6}\left(N_x-1\right)\left(2N_x-1\right)-{\overline{x}}^2\end{array}$$E30

The standard error ${\widehat{\sigma }}_x$ in the simulated parameter search estimates in Table 8 is of the same order of magnitude as that obtained from Eq. (30) with a sizeable discrepancy in the viscous friction coefficient which is possibly due to the limited B parameter sample space used. The initial temperature $T_i$ is determined from the resultant sample cost distribution illustrated in Figure 32, with details in Table 8, as

Figure 32.

Histogram of exploratory search costs.

$$T_i=\mathit{k\sigma }=1.654.$$E31

The anneal temperature scaling factor k is chosen [60] with a typical value of 10 on the basis of acceptance of parameter configuration costs, assumed normally distributed, worse than the present value by $3\sigma$ with a successful Boltzmann jump probability of

$$P=exp\left(-\mathrm{\Delta }E/T_i\right)=exp\left(-3\sigma /T_i\right)\approx 0.75$$E32

such that

$$k=-3/lnP=10$$E33

In this statistical gathering phase the best estimate ${\widehat{\mathbf{X}}}_{\mathit{\text{iopt}}}$ of the optimal parameter vector with cost ${\widehat{E}}_{\mathit{\text{iopt}}}$ is retained for initialization of the ensuing FSD algorithm if required. The cooling schedule is the most critical feature of the FSD method to guarantee global convergence and avoid trapping at local minima due to premature quenching. The update temperature algorithm consists of a sequence of temperature decrements and a condition to secure thermal equilibrium at each stage so that the average cost decreases in a uniform manner overall. The sequential temperature decrease

$$T_{k+1}={\alpha }_k\ast T_k$$E34

is based on the step reduction factor

$${\alpha }_k=exp\left(-\lambda \frac{T_k}{\sigma }\right)$$E35

with a typical value for λ of 0.8 and has been reported to work well in practice [24, 60] for a diversity of applications. The following numerical constants used in the FSD parameter extraction process have been obtained from heuristics and are known to give good performance of the algorithm with accurate estimates of global minimum convergence [60]. In practical applications the reduction process is typically lower bounded by

$$\mathit{min}{\alpha }_k=0.5$$E36

to prevent too rapid a reduction. Acceptance of large cost function jumps is a feature of the high temperature phase as in simulated annealing obviating the need for the ineffective gradient term in the iteration process and thus only a random search of parameter space is made initially for the first 10 external loops. At each temperature step the number of iterations performed is predefined at 15 times the dimensionality of parameter space (NDIM) to ensure equilibrium of the method. As the control temperature is reduced the volume of random search space is pruned, which is an attribute of the $\sqrt{2T}$ term of the FSD equation in Eq. (17), in accordance with a temperature dependent multiplier proportional to

$${\left(\frac{T_{k+1}}{T_i}\right)}^r,0.5\le \mathrm{r}\le 1.0.$$E37

This enhances convergence towards the global parameter optimizer X_opt and curtails non profitable random moves. Another useful feature of the FSD random search [24] is the application of a long tailed Lorentzian parameter distribution, which permits the occurrence of large cost jumps from local minima at low temperatures and consequent metastable state ${\widehat{\mathbf{X}}}_{\mathit{opt}}^i$ trap avoidance. At the termination stage of the FSD algorithm, when there is little observed change in the cost function over successive temperature steps, a controlled reheat phase is introduced with gradual temperature increase to reduce the risk of premature quenching and local minima trapping. The final phase of temperature decrease to the freeze condition results in convergence to the best approximation of the global minimiser.

9. Results obtained from FSD parameter identification

The extraction of the dynamical parameters for a typical brushless motor drive system with zero shaft load (NSL) inertia using the FSD method [26] is qualitatively illustrated in Figures 33 and 34 with returned estimates summarized in Table 10. The application of a uniform random exploratory search of parameter space initially is rewarded by the provision of a good estimate of the optimal parameter set for initialization of the FSD algorithm. The returned estimate ${\widehat{\mathbf{X}}}_{\mathit{\text{iopt}}}$ in Table 10, which has a cost ${\widehat{E}}_{\mathit{\text{iopt}}}$ equal to 15% of the exploratory sample mean $\overline{E}$, in this case is very close to the best possible estimate ${\overline{\mathbf{X}}}_{\mathit{opt}}$ of the optimal vector with the minimal quadratic penalty ${\overline{E}}_{\mathit{opt}}$ for the parameter tolerance range chosen. The best estimate ${\overline{\mathbf{X}}}_{\mathit{opt}}$ for the optimal parameter set can be obtained by visual inspection from the experimental I_fa cost surface in Figure 35 employing worst case quantized parameter step sizes in Table 8. This estimate ${\overline{\mathbf{X}}}_{\mathit{opt}}$, which will be referred to from here on as the ‘optimal’ vector for convenience and brevity of expression, is used as a reference against which the accuracy of the returned FSD parameter estimates can be judged. The proximity of ${\widehat{\mathbf{X}}}_{\mathit{\text{iopt}}}$ to ${\overline{\mathbf{X}}}_{\mathit{opt}}$, as determined from the error analysis in Table 11, is indicated by the degree of quantized resolution of parameter space with error differential

Figure 33.

FSD cooling sequence record.

Figure 34.

Iterative cost reduction sequence.

Figure 35.

MSE cost surface without shaft load.

Optimal Estimates obtained via Experimental Cost Surface in Figure 35
Quantized Parameter FC Response Surface Simulation Mesh Size $\left[2\delta J^L,\delta B^L\right]$
${\overline{J}}_{\mathit{opt}}$ =3.0797x10^−4 kg⋅m^2	${\overline{B}}_{\mathit{opt}}$=1.9207x10^−3 Nm/rad/sec	${\overline{E}}_{\mathit{opt}}$=5.0051x10^−2
Initial Exploratory Search and Returned FSD Estimates - Table 8
Sample Statistics	Ni = 200	σ = 1.654x10^−1	mean $\overline{E}$ =0.37
${\widehat{J}}_{\mathit{\text{iopt}}}$= 3.0638 kg⋅cm^2; =99.5%${\overline{J}}_{\mathit{opt}}$	${\widehat{B}}_{\mathit{\text{iopt}}}$ = 1.953x10^−3 Nm/rad/sec; =102%${\overline{B}}_{\mathit{opt}}$	${\widehat{E}}_{\mathit{\text{iopt}}}$ = 5.473x10^−2; =109%${\overline{E}}_{\mathit{opt}}$
Returned FSD Optimal Parameter Estimates
FSD Initialization: Xi	Ji = 82%Jm	Bi = 109%Bm
Observations
After First Temperature Step T1:   →		Main Lobe Capture as in Figure 36; J^(1)opt = 3.0341 kg⋅cm^2 = 98.5%${\overline{J}}_{\mathit{opt}}$; B^(1)opt = 2.199 mNm/rad/sec =114.5%${\overline{B}}_{\mathit{opt}}$; E^(1)opt = 7.346x10^−2 < σ
• After Fifth Temperature Step T5:  →; • First Reheat Cycle Completed:   →; • Total Number of Temperature Steps: →		Cost Reduction Ceases; 4 Temperature Steps {NF = 4}; J^(2)opt = 3.07058 kg⋅cm^2 = 99.7%Jopt; B^(2)opt = 1.9586 mNm/rad/sec =102%Bopt; E^(2)opt = 5.295x10^−2 < E^(1)opt < σ; 18 {First 10 used for Random Search}
FSD Global Convergence Estimate: ${\widehat{X}}_{\mathit{opt}}=\left\{J_{\mathit{opt}}^{\left(2\right)},B_{\mathit{opt}}^{\left(2\right)}\right\}$
Computational Details of FSD Algorithm in Figures 28 and 29
No. of External Temperature Loops with Random Search only mj = 10
No. of Parameter Adjustments per Temp. Step to reach Equilibrium uj = 15⋅NDIM = 30
No. of MSE Evals.→	Random: 405	Gradient: 387
Av. No. of BLMD Simul^ns per Temp. Step: →		(405+387)/18 = 44
Total Simulation Time for FSD Algorithm: →		20546 secs for a 486-DX-66 MHz CPU
Execution Time per BLMD Simul. Trial: →		20546/(405+387) ≈ 26 sec

Table 10.

Summary of returned FSD parameter estimates.

BLMD Parameter	Motor Inertia Jopt; $\left(\frac{{\widehat{\mathrm{J}}}_{\mathrm{opt}}-{\overline{\mathrm{J}}}_{\mathrm{opt}}}{{\overline{\mathrm{J}}}_{\mathrm{opt}}}\right)\times 100\%$	Damping Factor Bopt; $\left(\frac{{\widehat{\mathrm{B}}}_{\mathrm{opt}}-{\overline{\mathrm{B}}}_{\mathrm{opt}}}{{\overline{\mathrm{B}}}_{\mathrm{opt}}}\right)\times 100\%$	Global Cost Eopt; $\left(\frac{{\widehat{\mathrm{E}}}_{\mathrm{opt}}-{\overline{\mathrm{E}}}_{\mathrm{opt}}}{{\overline{\mathrm{E}}}_{\mathrm{opt}}}\right)\times 100\%$
Exploratory Search	−0.515%	1.682%	9.348%
FSD Method; Post Temp Step T1; Post Temp Step T5	−1.479%; −0.296%	14.49%; 1.976%	46.77%; 5.808%
Coherence of BLMD Simulation Waveforms with Experimental Test Data
NSL - Correlation; Coefficient ρ	Current Feedback Ifa; 94.62%	Current command Ida; 95.815%	Cur. Controller o/p Vca; 94.9%

Table 11.

Error Analysis for Returned FSD Parameter Estimates.

Shaft Velocity Correlation Coefficient V_ωr = 97.9%.

$$d{\widehat{\mathit{X}}}_{\mathit{\text{iopt}}}={\left[2\delta J^L,\delta B^L\right]}^{\mathrm{T}}.$$E38

The effectiveness of the FSD method in achieving parameter ‘optimality’ over the FC undulating cost surface is demonstrated by deliberately initializing the search far from the global minimum in the neighborhood of a local minimum at Xi as given in Table 10.

A scatter diagram of iterative search costs levied by the FSD process is portrayed in Figure 36 and contrasted with a cross section of the FC cost surface shown in Figure 37 for variable J and fixed $B={\overline{B}}_{\mathit{opt}}$. Application of a Boltzmann probability step transition in surmounting significant cost barriers, with cost elevations of 2σ above $\overline{\mathrm{E}}$ which conceal the global minimum, is evident in the iterative search improvement process for a global optimum. Thus a means of escape from local minimum capture, which causes problems for traditional identification methods with non-optimal convergence, is provided for parameter extraction in an effort to secure the least mean squares (LMS) estimates.

Figure 36.

FSD parameter search with local minimum escape.

Figure 37.

Acquisition of Global Minimum from experimental I_fa cost surface.

In the early stages of the cooling schedule the temperature reduction is matched by a rapid reduction in iterative cost which eventually saturates as the parameter estimates approach the optimal LMS values. After the first temperature step the estimation cost, which is less than the exploratory $\sigma$, is captured within the main lobe of the response surface containing the global minimum. Subsequent progress in cost improvement is minimal which is mainly due to the fact that the parameter estimate at X⁽¹⁾_opt is relatively close to the global minimizer and that random searching is employed for the first ten temperature steps without the assistance of a gradient search. In addition to this the selectivity of the cost surface to B parameter variation is poor along the valley floor, which is flat near the global minimum as illustrated in Figure 35 as per [13, 28], making it difficult for global minimum convergence.

During this phase of cost immobility a reheat cycle is introduced [24], as part of the FSD cooling schedule shown in Figure 29 to overcome local minimum trapping, which is irrelevant in this case as the current estimate is circumjacent the optimal value. After the fifth iteration further cost reduction ceases with the global convergence estimate given by

$${\widehat{\mathbf{X}}}_{\mathit{opt}}={\mathbf{X}}_{\mathit{opt}}^{\left(2\right)}.$$E39

A second reheat phase followed by freeze conditions does not improve the parameter estimate. The total execution time on a 486-DX-66 MHz processor for the FSD method in this parameter identification procedure amounted to 20546 seconds with a computation burden of 792 cost function evaluations over 18 temperature steps. This translates on average into 44 BLMD model simulation runs per temperature step with a processor runtime of 26 seconds per simulation trial. This puts into perspective the computational intensity and the cost in terms of CPU time of this optimization technique. A breakdown of the overall search effort deployed in the FSD process projects into 405 random search calculations with a punitive computational overhead of 387 evaluations necessary for anticipated downhill movement factored in only for the last 9 temperature steps. This overhead is manifested, after the random search epoch during which the parameter estimates have reached a minimum potential, in the lengthy reheat and condensation phases shown in Figure 33. The inadequacy of the gradient search procedure in this instance, with approximate global convergence already achieved, is compounded by the tendency of the FSD method to dither about in a region of parameter space where the response surface has poor selectivity in the B parameter with little iterative progress in cost reduction expected. The comparison of BLMD step response simulation with experimental drive test data, at internal observation nodes based on returned FSD optimal estimates, provides an excellent fit in terms of frequency and phase coherence as per Figures 38–40 pertaining to BLMD current control.

Figure 38.

ZSL FC Simulⁿ via MFSD estimates.

Figure 39.

ZSL- current demand Simulⁿ.

Figure 40.

ZSL-current controller Simulⁿ.

These simulated step responses are almost identical to those based on returned estimates for the modified form of the FSD method (MFSD), with zero inertial load conditions (NSL), discussed in Section 11 below. These simulation traces are displayed together along with other waveforms for MFSD comparison purposes in Figures 38–41 and Figures 42–49 at different shaft inertial load bearing conditions. The goodness-of-fit measure is determined by the degree of correlation, expressed as the correlation coefficient in Table 11, which exists between the simulated and experimental test data. This value is almost 95% for all three internal node observations which confirms the success of the FSD method as an accurate method of parameter extraction besides providing BLMD model validation. Further confidence in the FSD method is assured by the degree of correlation in Table 11 of BLMD shaft velocity simulation with test data as shown in Figure 41.

Figure 41.

ZSL – MFSD Shaft Velocity Simulation.

Figure 42.

FC simulation via MFSD estimates.

Figure 43.

MSL - current demand simulation.

Figure 44.

Current controller o/p simulation.

Figure 45.

FC–LSL Simulation via MFSD Estimates.

Figure 46.

LSL - current demand simulation.

Figure 47.

LSL - current controller simulation.

Figure 48.

MFSD – MSL Shaft Velocity Simulation.

Figure 49.

MFSD – LSL - Shaft Velocity Simulation.

10. Improved FSD method for motor parameter identification

A drawback with the general form of the FSD method described in Section 8 is the usage of a fixed number of external loops at high temperature for random searching only of parameter space without application of the gradient technique. The lack of adaptation of the method in its present form to the iteration statistics garnered during the initial high temperature phase gives rise to lengthy run time of the FSD algorithm. This results in an excessive number of cost calculations, which can be very expensive in central processing unit (CPU) time resources, accompanied by at least one reheat phase before termination of the FSD method without cognisance of the gains made at each temperature step and the shape of the cost surface based on motor feedback current.

It is evident from Figure 34 that significant reductions in the objection function are obtained during the first couple of high temperature steps only, with very little improvement thereafter. A more lucrative strategy [27] can be based on adaptation of the FSD method during the high temperature random search phase in the absence of good initialization of the parameter vector and on the initial exploratory search statistics where the maximum value E_max, mean $\overline{E}$ and standard deviation σ of the resultant cost distribution are known.

The initial exploratory search statistics, for three known cases of motor shaft inertia to be ‘identified’, are listed in Table 12. The histograms associated with the initial exploration of inertial parameter space, along with those for the related costs, are shown in Figures 50–53. These histograms possess the same uniform randomness attributes as that for the zero shaft load inertia case discussed in Section 8 with cost distribution depicted in Figure 32.

Total Nominal Inertial; Shaft Load; Jnom ($\mathrm{kg}.{\mathrm{cm}}^2$)	No Shaft Load; NSL - (0.0); 3.0	Medium Shaft Load; MSL - (9.304); 12.304	Large Shaft Load; LSL - (17.822); 20.822
Tolerance: ±ΔXnom = [±20% Jnom, ±10% Bnom]		Sample Size: $N_i=200$
Quantized Parameter Space: Sample Percentage	9.542%	2.315%	1.19%
Max Cost ${\mathrm{E}}_{\max }$	0.694	0.655	0.653
Mean Cost $\overline{E}$	0.367	0.444	0.457
Min Cost ${\widehat{E}}_{\mathit{\text{iopt}}}$	$5.473\times {10}^{‐2}$	$2.054\times {10}^{‐2}$	$1.497\times {10}^{‐2}$
Standard Deviation σ	$1.654\times {10}^{-1}$	$1.651\times {10}^{-1}$	$1.694\times {10}^{-1}$
Initial Optimal Estimates; ${\widehat{J}}_{\mathit{\text{iopt}}}\left(\mathrm{kg}.{\mathrm{cm}}^2\right)$; ${\widehat{B}}_{\mathit{\text{iopt}}}\left(\mathrm{Nm}.{\mathrm{rad}}^{‐1}.sec\right)$	3.064; $1.953\times {10}^{‐3}$	12.121; $1.953\times {10}^{‐3}$	20.767; $2.033\times {10}^{‐3}$

Table 12.

Statistical estimates for initial exploratory search of parameter space.

Figure 50.

Histogram of MSL-J search values.

Figure 51.

Histogram of MSL exploratory costs.

Figure 52.

Histogram of LSL-J Search Values.

Figure 53.

Histogram of LSL exploratory costs.

It is evident that all the tabulated exploration costs fall within $\pm 3\sigma$ of the sample mean $\overline{E}$ with worst case estimates E_max less than $\overline{E}+2\sigma$, for different shaft load inertia as illustrated in Figure 37, which according to Chebyshev’s theorem of measurements [13, 28, 61] should account for at least 89% of all costs. Furthermore exploratory parameter estimates ${\widehat{\mathbf{X}}}_{\mathit{\text{iopt}}}$ with costs less than the rms deviation σ about the mean are very close to being optimal and are in the ambit of global minimum convergence given the shape of the FC response surface. Such estimation costs, which are less than $\overline{E}-\sigma$, are confined to the main lobe of the objective function which has for an apex the global minimizer ${\overline{\mathbf{X}}}_{\mathit{opt}}$. For the purpose of validation and demonstration of the effectiveness of the MFSD method as a parameter extraction tool in identifying various motor shaft inertial loads the exploratory least squares parameter estimates in Table 12 are ignored in the initialization process.

Instead the method is deliberately initiated far from the global estimates in all three cases at

$$X_{\mathit{init}}={\left[82\%J_{\mathit{nom}},109\%B_m\right]}^{\mathrm{T}},$$E40

based on parameter tolerance bounds given in Table 12 with initial temperature given by

$$T_i=T_{\mathit{\text{init}}}=10\mathrm{\sigma },$$E41

where the objective function value is large and in the vicinity of a local minimum. The best optimal parameter estimates are obtained separately by inspection from quantized parameter FC cost surface simulations, utilizing the same experimental FC target data, for comparison purposes with the MFSD method as shown in Figures 35,54 and 55 for different inertial loads. The modified version of the FSD algorithm uses a discretized parameter space with step sizes as listed in [13] and incorporates the following improvements:

1. proceed with normal random probing of parameter space for the first high temperature step as in Figures 28 and 29 with the best parameter estimate retained. At the end of this iterative sequence a ‘greedy’ search is performed, from the current best estimate obtained thus far, along each of the parameter directions in turn using the gradient method with a quantized parameter step size to take advantage of any possible further cost improvement that may accrue. The resultant optimal parameter estimates are trimmed to the nearest step size.

2. If the current Least Mean Squares (LMS) error is less than $\sigma$ then the alternate application of hill descent with random search is pursued during subsequent temperature steps in the iterative improvement process towards optimality. This approach is adopted because the best estimate available with cost σ most likely resides within the capture zone of the main cost lobe as depicted in Figure 37 for example. If the converse is true then random searching is continued as in step (a) until sufficient cost reduction has been achieved.

3. A gradient search at the termination of each temperature step is maintained for further gain in cost reduction. This has the effect of forcing and hastening the convergence of the parameter estimates towards optimality. If however there is no cost improvement for several temperature steps, typically 3 and pending step (b), the gradient search is ceased and the FSD process enters the termination stage after the fourth step which corresponds with the onset of the freeze condition in the normal FSD method. The termination process can be speeded up by halting the algorithm at this point without the necessity of the freeze condition as subsequent parameter convergence information bears out in Table 13. This three step stalled cost reduction phase allows for improved convergence in the B parameter estimate along the valley of the response surface, where it is very flat in this parameter direction, without undue expenditure of computational effort. Secondly the anneal temperature rises during this interregnum, in line with the reheat phase of the FSD method, nullifying the possibility of false minima trapping of parameter estimates.

Figure 54.

MSE surface with medium shaft load.

Figure 55.

MSE surface with large shaft load.

Total Nominal Shaft Load Inertia; Jnom ($\mathrm{kg}.{\mathrm{cm}}^2$)	Figure 35 - NSL; No Shaft Load; Jm = 3.0	Figure 54 - MSL; Medium Inertia; 12.304	Figure 55 - LSL; Large Inertia; 20.822
Parameter Initialization: [Ji = 0.82Jnom, Bi = 1.09Bm] with Bm = 2.14$\times {10}^{‐3}$
Initial Cost Ei	0.414	0.590	0.616
Returned FSD Optimal Parameter Estimates
${\widehat{\mathrm{J}}}_{\mathrm{opt}}$ ($\mathrm{kg}.{\mathrm{cm}}^2$); ${\widehat{\mathrm{B}}}_{\mathrm{opt}}\left(\mathrm{Nm}.{\mathrm{rad}}^{‐1}.sec\right)$	3.07965; 1.95863$\times {10}^{‐3}$	12.13966; 2.00082$\times {10}^{‐3}$	20.8375; 1.95863$\times {10}^{‐3}$
Min Cost ${\widehat{E}}_{\mathit{opt}}$	5.21677 × 10^−2	2.20282 × 10^−2	1.30295 × 10^−2
No. of Temp Steps to reach ${\widehat{E}}_{\mathit{opt}}$	2	2	2
No. of Func_Evals; Random Search; Gradient Search	110; 45; 65	109; 45; 64	103; 45; 58
Total Time (sec)	2834	5562	6499
Average Iteration Time tITER	1417	2781	3250
Simul^n Time/Func. Eval. tSIM	25.76 s	51.03 s	63.10 s

Table 13.

Modified FSD method of motor shaft inertial load parameter extraction.

Once the motor friction coefficient B has been resolved further pruning of CPU run time can availed of if only the inertial parameter is to be identified, for which the FC response surface is more sensitive near the global minimum, by abrupt termination of the FSD method after two non profitable temperature steps.

11. Parameter convergence results for the modified FSD method

The optimized parameter estimates of the motor dynamics, returned by the modified FSD method [25], with convergence details for three cases of shaft load inertia to be identified are summarized in Table 13. The number of iterative temperature steps and functional evaluations required, for the motor parameter identification process to reach minimum potential, in each of the three cases of shaft load inertia are almost identical for the modified FSD method with considerable savings effected in computational effort over the conventional FSD approach.

The anneal temperature profiles realized during the MFSD extraction process of the three shaft inertial loads are depicted in Figures 56–58. Global optimality, to within the limits of the quantized parameter step sizes listed in [13], is achieved in all three cases in at most two temperature steps as shown in Figures 59–61. The CPU runtime however increases with the mechanical time constant τ_m given in [9], for a similar number of motor simulation trials as indicated in Figure 62, in symphony with the time duration of the observed FC target data used as an argument in the objective function formulation.

Figure 56.

NSL cooling temperature history.

Figure 57.

MSL cooling temperature history.

Figure 58.

MSL cooling temp. Sequence.

Figure 59.

Iterative NSL cost reduction.

Figure 60.

Iterative MSL cost reduction.

Figure 61.

Iterative LSL cost reduction.

Figure 62.

MFSD computation time.

While in general increasing the number of FC cycles improves the accuracy of the extracted global parameter estimates, measurement constraints limited the size of the FC transient data record to 10 cycles used in the MSE cost formulation. The length of each data training record is fixed at 4095 sample points (N_d) with a normalized time duration of approximately 10 machine cycles for reference purposes and response surface comparison with details given in Table 3.

Response surface simulation, although computationally expensive, provides an alternative route of accurately obtaining the optimal parameter vector by means of inspection of the surface minimum cost. It can thus be used as a yardstick by which the overall convergence performance of the MFSD method can be contrasted over a range of motor shaft load inertia. The cost simulations are based on an initial crude parameter mesh size given in Table 3 with a 1 μs time step Δt in all cases as shown in Figures 35,54 and 55 with further refinement down to quantized step sizes necessary in the vicinity of the global minimum for resolution of the optimal parameter set.

The optimal shaft inertia values extracted by this approach are used in the error analysis given in Table 5 as a reference by which the accuracy of the FSD optimal estimates are gauged. The parameter estimates extracted by the FSD optimization technique are very accurate in the shaft inertia only with fractional percentage relative errors achieved in global convergence performance as indicated in Table 5.

The errors appear to decrease with increasing inertial load coupled with a general trend in the global cost reduction. The error performance in B parameter estimation is ‘poor’ by comparison and is responsible for the increase in relative global cost in column 2 of the error summary. The wider error margins can be inferred from the deployment of the larger quantized relative step size listed in [13] and attributed to poor selectivity of the cost surface, in the friction coefficient, in the region of the global minimum as discussed in [13]. The quality of the MFSD probe of parameter space in obtaining the global estimates and evading local minimum capture is demonstrated in Figures 63–65 by the asterisked search costs over cross sections of the simulated FC response surfaces for contrast.

Figure 63.

FSD identification of rotor dynamics.

Figure 64.

FSD identification of inertial MSL.

Figure 65.

FSD identification of inertial LSL.

The accuracy of the returned MFSD estimates can also be checked by employing these optimal parameter vectors in BLMD dynamical simulation and resultant comparison with experimental motor drive step response data. The model simulation traces are compared with observed test data, on the basis of accuracy of fit at critical internal nodes of the drive system, for the dual purpose of validation of the MFSD identification method and enhancement of BLMD model confidence which are essential intrinsic component features of system identification [35].

The goodness of fit of the simulation traces, at various BLMD model observation nodes in [13] in terms of frequency and phase coherence, with sampled test data is obvious from Figures 38–40 and from Figures 42–47 for BLMD feedback current, current demand and compensator outputs at different inertial loads. This measure of trace coherence, indicative of MFSD parameter extraction accuracy, can be gauged by the excellent correlation coefficient for the fixed amplitude swept frequency waveforms listed in Table 14 (Figures 38–41,42–49). The accuracy of the MFSD method is further substantiated in Figures 41,48 and 49 by the correlation of the model shaft velocity step response with BLMD output test data.

Nominal Shaft Inertial Load; Jnom ($\mathrm{kg}.{\mathrm{cm}}^2$)	No Shaft Load (NSL); Jm = 3.0	Medium Inertia; (MSL); 12.304	Large Inertial Load; (LSL); 20.822
Current Feedback; Ifa	Figure 38; 94.6%	Figure 42; 98.6%	Figure 45; 99.4%
Current command; Ida	Figure 39; 95.8%	Figure 43; 97.8%	Figure 46; 99.7%
Current Controller o/p; Vca	Figure 40; 94.9%	Figure 44; 99%	Figure 47; 98.7
Motor Shaft Velocity; Vωr	Figure 41; 97.9%	Figure 48; 99%	Figure 49; 93.8%

Table 14.

BLMD simulation trace coherence correlation coefficient ρ in [13].

12. FC response surface selectivity with MFSD application

The above application of the modified version of the FSD method in inertial mass J identification evolved from the nature of the FC objective function deployed. The fixed sample number N_i (200) of random searches during the exploratory phase, or the first temperature step employing 15*NDIM random searches with user supplied initialization only and an absent exploratory phase, tended to give less coverage of quantized parameter space for the cases with increased shaft inertia presented in Table 12. It therefore seems more difficult to find a good starting vector, as it takes longer to probe a sufficient volume of parameter space, from which to anchor all subsequent searches. The modified method would, it appears, then terminate prematurely if trapped in a false minimum potential well after four anneal steps without adequate sampling of the parameter domain thus generating uncertainty as to the quality of returned estimates. However this is not the case when the topography of the cost surface in the observed FC variable is examined for different shaft inertial loads.

The number of feedback current cycles is fixed, at about ten in Table 3, as this influences the FC response surface “nuclear capture cross section” [62] shown in Figure 66. This can be defined as the waist of the main lobe containing the global minimum ‘energy’ E_opt, with value less than the exploratory search σ, and the penalty costs of all excursions into the parameter space kernel forming a nucleus about X_opt with values less than the next excited level E₁ or local minimum above the ground state E_opt. The probability of capture P_c of a random search can be formulated as the ratio of the capture cross section to the tolerance band of parameter variation in [13] about the nominal value x_m as

Figure 66.

FC cost surface ‘nuclear’ capture zone.

$$P_C=\frac{\Delta x_{\mathrm{cap}‐\mathrm{x}‐\text{sect}}}{x_{\max }-x_{\min }}=\frac{td{x}^L}{2md{x}^L}=\frac{t}{2m}$$E42

where it is assumed that the least important B parameter value has been resolved. Thus for a uniform random search the impingement of target parameter space, given by NP_c, increases directly with the number of trials N. Once the barrier potential E₁(X) to the kernel of capture is breached at any stage during the MFSD iterative improvement process, the related parameter vector is then installed as the best estimate by which all subsequent searches are adjudicated on for any further cost decrease. At the end of an anneal step a gradient search is performed which forces the best estimate to date, if trapped in the main lobe region, to converge towards the global minimum.

The capture cross section percentages, expressed in terms of allowed parameter tolerance and quantization step size, for the three cases of motor shaft inertia are presented in Table 15 along with details of search cost capture during implementation of the modified FSD method. The tabulated cross section percentages stabilize at about 30% for medium to large mass loading which shows that the probability of reaching the target zone is independent of the number of searches conducted owing to the shape of the optimal cost surface sections employed. This removes any uncertainty as to the effectiveness of the FSD method of entering the global region and possible lockup in an excited state in the event of a short anneal sequence where a large quantized parameter space is concerned.

Capture Cross Section; at Bopt	No Shaft Load; Figure 35; 38.97%	Medium Inertia; Figure 54; 30.15%	Large Inertia; Figure 55; 29.69%
FSD Target Acquisition
Exploration Phase; Trapped States %	127 out of 200; 63.5%	85 out of 200; 42.5%	71 out of 200; 35.5%
First Anneal Step; Random Search Only: %Trapping	15 ex 30; 50%	16 ex 30; 53.33%	12 ex 30; 40%
First 2 Temp Steps; Random + Mixture: %Capture	81 out of 110; 73.64%	75 out of 109; 68.81%	78 out of 103; 75.73%

Table 15.

Estimation of % capture cross section from simulated FC cost surfaces.

Also the various FSD entrapment estimates at different phases of the search algorithm provide a rough independent confirmatory measure of the probability of capture where a uniform random search of parameter space is conducted. The percentage of all search excursions into the optimal parameter nucleus, in both the optional exploration phase and the mandatory random initial anneal step, exceed the cross section estimates in all cases and thus enhances user confidence in the modified method. The rationale behind this approach is based on the FC cost surface topography, which has a ravine like capture region shown in [13], with both J and B parameters considered. This target zone, which has essentially parallel contour lines in the direction of the ‘line minimum’ in [13], can be condensed into a single dimension in the percentage ratio carve up of parameter space for target acquisition calculations in the J parameter as shown Figure 66.

Further confidence in the capture cross section estimate for different shaft inertial loads can be gained by employing fixed length data records in the observed current feedback I_fa. This is based on ten FC machine cycles as the target reference for FC cost surface simulation in the J parameter at B_opt for estimation purposes as summarized in Table 16 with data sampling rates f_s and displayed in Figure 67.

Figure 67.

C_X_S variation with inertia.

The probability of capture can be more accurately defined by the weighted contribution of all tabulated cross section estimates at different inertial loads as

$$P_C=\frac{{\sum }_i{J}_i\ast \mathrm{C}\_\mathrm{X}\_\mathrm{S}}{{\sum }_i{J}_i}=0.31.$$E43

The adoption of a fixed threshold step size δJ^L during MFSD parameter extraction results in improved parameter resolution, given by the reduction in the variability estimate as

$$V_x^L=\frac{\delta x^L}{x},$$E44

with increased shaft inertial loads as shown in Table 9.

This percentage reduction results in better accuracy of the returned J parameter values for increased shaft inertial loads with smaller relative errors as indicated by error analysis in Table 5. The selectivity of the FC objective function in parameter identification can be improved somewhat by increasing the length of the I_fa sampled data record used as the target reference in the least squares cost formulation. The convergence domain is reduced with the data record length in an inverse time relationship thus improving the accuracy of the returned parameter estimates. This is articulated by the accompanying reduction in capture cross section thus restricting the search for optimality to a smaller kernel within the parameter nucleus surrounding X_opt. The raison d’être is the eclipsing of the finite length motor speed exponential ramp up transient by the ever lengthening presence of steady state speed saturation. The capture cross section variation for different simulated data file lengths, with a fixed sampling rate of 50 kHz and ± 20% parameter tolerance, is illustrated in Table 17 for nominal rotor inertia J_m anchored at ${\overline{B}}_{\mathit{opt}}$.

Ifa [fs]	50 kHz	37 kHz	31 kHz	27 kHz	24 kHz	21 kHz	17 kHz	15 kHz
Inertia; $\mathrm{kg}.{\mathrm{cm}}^2$	3.0	6.0	9.0	12.0	15.0	20.0	30.0	40.0
C_X_S	36.3%	33.9%	32.3%	31.7%	31.2%	30.5%	30.2%	30.3%

Table 16.

% capture cross section (C_X_S) estimates from I_fa target data.

The accompanying graph in Figure 68 shows that the capture cross section of the optimal parameter kernel decreases with increasing data block size and ultimately stabilizes to a constant value for modest lengths which is indicative of the fact that the motor has reached maximum shaft speed. The increased selectivity is accompanied by a plurality of excited states, shown in Figure 69, which eventually level off in number as the shaft speed reaches its maximum value. This manifestation is due to the contributory interference effect of the two frequency modulated sinusoids, during the machine speed step response, in the penalty function construct as explained in [13]. The selectivity is obtained by fitting a quadratic polynomial, expressed in terms of the parameter J, to the main lobe of the cost surface cross section centred on the vertex ${\overline{\mathbf{X}}}_{\mathit{opt}}$ as shown in Figure 69 with

Figure 68.

C_X_S versus data length.

Figure 69.

Local minimum proliferation.

$$E\left(J\right)-{\overline{E}}_{\mathit{opt}}=b_0{\left(J-{\overline{J}}_{\mathit{opt}}\right)}^2$$E45

and least squares coefficient b₀ determined [9, 10]. The selectivity measure is based on the first order variation δE(J) of the penalty cost function about the global minimum estimate ${\overline{X}}_{\mathit{opt}}$ in terms of the worst case parameter step size δJ^L as

$${\left.S_J^E\right|}_{{\overline{X}}_{\mathit{opt}}}=\frac{\mathit{\delta E}}{\mathit{\delta J}}=b_0\delta J^L.$$E46

This measure indicates a quasi linear dependency with increased data record length N_d as detailed in Table 9 and shown in Figure 70. Further increase in data record length causes the capture cross section to tend asymptotically to zero, with the nearest local minimum to the global entity approaching the same minimum energy as the ground state E_opt, due to the gradual onset of steady state conditions. The increased data length results in a proliferation and clustering of genuine local minima, other than noise induced cost surface granularity, about the global extremum as the initial transient step response is swamped by the onset of steady state motor speed.

Parameter Vector Fulcrum: $X={\left[J_m,{\overline{B}}_{\mathit{opt}}\right]}^T={\left[3.0\times {10}^{-4},1.921\times {10}^{-3}\right]}^T$
No. Samples	4120	6280	9000	12000	24000
Shaft Speed RPM	1953	2441	2686	2808	2930
C_X_S	36.3%	19.69%	12.56%	9.24%	5.1%

Table 17.

Response surface selectivity improvement with target data length.

Figure 70.

FC cost surface selectivity.

It is noteworthy that a judicious application of the Fast Fourier Transform (FFT) to the FC target data, after the initial transient has partially expired, provides an estimate of the motor speed variation with residual data record length as shown in Figure 27. The BLMD shaft speed can be approximated by the FFT of the FC step response over a short time span at the end of each FC data sequence where machine speed is quasi constant and shaft velocity fluctuations have almost disappeared. The exponential buildup of the motor shaft speed with data record length/time in Figure 27 is indicative of the characteristic response to the torque demand step input [9, 10].

13. Conclusions

The FSD optimization technique has been shown, although computationally intensive, to be a very accurate and effective parameter identification method over a noisy cost surface with embedded local minima. This novel method avoids the convergence difficulties associated with the application of classical optimization techniques during parameter extraction by providing an anneal temperature related Boltzmann probability of escape from local minimum capture. The accuracy of the FSD returned parameter estimates is independently confirmed with values found to be in substantial agreement with those extracted through cost surface simulation, with worst case relative error in the damping parameter B of 2%, for known BLMD dynamics. Further confidence enhancement in the FSD identification strategy is provided by the degree of correlation (∼94%) of BLMD model simulations, based on returned FSD parameter estimates, with observed test data. The application of the Powell Conjugate Direction method, as an alternative competitive identification strategy to the FSD method, with a simple unimodal velocity cost function has been discussed. The identified BLMD dynamics, though reasonable by comparison with known nominal values in that relative errors are less than 6% in Table 4, are not as accurate as those returned by the FSD method. This discrepancy is apparent for BLMD waveform simulation over long spans, using returned PCD parameter estimates for three known cases of inertial load to be identified, with poor correlation estimates spanning 20–70% as per Table 6 for FM related current traces. The reasons, relating to cost surface selectivity and returned parameter accuracy, underpinning the choice of target data in the MSE objective function have been discussed. The obscuration of shaft velocity transient response detail by the onset of steady state conditions, as a consequence of lengthy data records, results in impairment of the PCD method in terms of reduced cost surface selectivity and loss of returned parameter accuracy with increased convergence metric in Table 7. The reverse trend emerges during FSD optimization, with better selectivity and improved parameter accuracy with smaller capture cross section, for current feedback usage with this being the preferred target data choice in MSE cost formulation.

The effectiveness of the modified version of the FSD algorithm has been demonstrated in the extraction of the BLMD dynamics with enhanced speed in global convergence which is less than 15% of the runtime for the original FSD algorithm based on zero inertial shaft load conditions. The accuracy of the method is verified with initialization far from the global minimum at the edge of the parameter tolerance band and supported by the returned results with relative errors less than 0.3% as per Table 5 for three known cases of motor shaft inertia. The correlation accuracy of BLMD waveform simulation with measured data is greater than 94% in Table 14 for all cases of returned inertial parameter estimates by the FSD method. The computation intensity of the method can be reduced by adopting coarser quantization step sizes, instead of the specified noise threshold value, within an acceptable bounded parameter relative error. This improved FSD method can be potentially used, as an effective optimizer in the dynamical parameter extraction phase, to facilitate autotuning during embedded system commissioning.

Acknowledgments

The author wishes to acknowledge

1. Eolas (Science Fundation Ireland) – The Irish Science and Technology Agency – for research funding.

2. Moog Ireland Ltd. for brushless motor drive equipment for research purposes.