Open access

Mathematical Analysis for Response Surface Parameter Identification of Motor Dynamics in Electric Vehicle Propulsion Control

Written By

Richard A. Guinee

Submitted: 23 August 2012 Published: 19 December 2012

DOI: 10.5772/54483

From the Edited Volume

New Generation of Electric Vehicles

Edited by Zoran Stevic

Chapter metrics overview

2,856 Chapter Downloads

View Full Metrics

1. Introduction

This chapter addresses the topographical examination of various mean squared error (MSE) cost surface structures and selecting the most suitable MSE fitness function for accurate brushless motor drive (BLMD) dynamical parameter system identification (SI) of BLMD shaft load inertia and viscous damping for electric vehicle controlled propulsion. The parameter extraction procedure employed here is in the offline mode for optimal drive tuning purposes during the installation and commissioning phase of embedded BLMD systems in high performance electric vehicle torque, speed and position control scenarios. Two types of penalty function, based on the transient step response of the permanent magnet (PM) motor shaft velocity and its stator winding current feedback in torque control mode [1,2], are examined here for arbitration of a suitable choice of cost objective function as the response surface in the accurate extraction of the BLMD dynamics. The choice of a particular MSE cost surface as an objective function in BLMD load parameter identification is motivated by the need for reliable tuning of the proportional and integral (PI) term settings during the drive installation phase for controller robustness and optimal performance in adjustable speed drive (ASD) or torque controlled embedded PM motor applications for electric vehicle propulsion. This chapter will focus on the mathematical analysis of embedded motor drive dynamical parameter identification over an MSE multiminima response surface with the following key results obtained:

  1. the development of a novel quadratic mathematical model approximation for the investigation of the (i) nature of the MSE objective function and (ii) existence of a bounded MSE global minimum stationary region, based on transient step response motor current feedback signals, for mechanical parameter identification in sensorless drive torque control of electric vehicles.

  2. the examination of the phenomenon of multiminima proliferation in the MSE cost formulation due to target data choice and ‘noisiness’ arising from evaluation of pulse width modulated (PWM) edge transition times during BLMD simulation [1,2].

  3. the measurement of cost surface selectivity based on shaft velocity and current feedback target data and the decision favouring the choice of the latter data training record as the target function for dynamical parameter identification

  4. the development of a novel parameter quantization metric to overcome cost surface ‘noisiness’, arising from computational uncertainty in the simulated PWM edge transitions [1], for avoidance of local minima trapping in the MSE cost surface during identification of the BLMD dynamics.

  5. the development of a novel parameter convergence radius measure of encirclement of the cost surface stationary region global minimum, arising from the parameter quantization metric in (d), for determination of the bounds of accuracy that can be imposed on the returned estimates of the global optimum dynamical parameter vector during BLMD identification.

1.1. Motivation

BLMD control tuning is necessary during the commissioning phase of embedded drives applications, for accurate torque and speed control in electric vehicle propulsion systems setup [3,4], accurate robotic end effector [5] or CNC tool positioning [6], where detailed apriori knowledge of expected drive load inertia and friction parameters are unknown to the electric motor drive supplier/manufacturer in the intended application beforehand. The choice of ASD [7] in high performance industrial applications, such as a small electric vehicle [4], robot manipulator [8, 9, 10] or machine tool feed drive [6,4], is usually based on consideration of a BLMD manufacturer’s catalogued specifications, relating to drive performance capabilities and limitations, by the customer or embedded drive equipment designer/manufacturer. The BLMD selection is often done independently of the motor drive manufacturer by the equipment designer for reasons of embedded systems design confidentiality and second sourcing of matching drive equivalents from different manufacturers for the purpose of cost reduction and embedded product protection from obsolescence via alternative drive substitution. The range of motor sizes available and spread of possible BLMD embedded applications has resulted in the provision of flexible drive tuning facilities with either manual or autotuning features [11] by motor drive manufacturers as a sales and marketing expedient to embedded equipment designers. This flexible approach to drive tuning policy eliminates the need for the BLMD manufacturer to participate in the detailed design of embedded drive applications except in the provision of motor drive systems with high output torque and speed ranges to cater for a range of anticipated high performance applications [12,8,13]. BLMD systems with high peak current capability and fast response times due to low PM rotor mass are designed [6, 7] to handle large inertial load torques [14] experienced in robotic applications [4, 5] and electric vehicle propulsion systems, with a no-load to full-load inertia variation [9] of 10 is to 1. It is in response to this background of applications diversity, regarding the particular design details of embedded drive products about the size of inertial loads and friction coefficients encountered [4,6,9], that the present work on cost surface analysis for parameter extraction in electric vehicle control is directed from a motor manufacturer’s perspective.

Since the possible variation in the load dynamics of an intended BLMD application is unknown at the outset the initial task here for an end user is to identify the actual load inertia and friction coefficients experienced during startup of a given embedded drive in the offline mode for robust PI controller tuning [15]. In this scenario the customer has the flexibility of manually tuning the BLMD speed loop, which is provided as a PI adjustment option along with procedural details for tuning by the motor manufacturer, during the setup and commissioning phase for a particular ASD application. The challenge then posed for the motor drive designer in this instance is the provision of an automated tuning facility for the velocity or torque loop during the commissioning stage thus eliminating the need for any manual input by the customer. This feature requires the identification of a fixed embedded load configuration during setup and subsequent automated optimal configuration of the velocity controller PI terms [15]. In the absence of embedded load information the cost surfaces and identification methods investigated here focused on inertial load spreads for vehicular and robotic applications [9] of up to ten times the inherent motor shaft inertia as recommended by the BLMD manufacturer for the drive [16] modelled in [1].

The concept of a simulated cost surface is developed here [17] as an objective function to facilitate parameter extraction of the installed drive dynamics, during offline BLMD system identification, with MSE minimization. This methodology provides useful insight into the nature and formulation of the most suitable MSE objective function to be minimized, based on actual drive experimental test data available and BLMD model simulation, coupled with an effective system identification (SI) strategy for accurate motor parameter extraction [18]. This approach can also be used as an alternative means of providing the optimal set of extracted parameter estimates from inspection of the global minimum location on the simulated cost surface with embedded local minima. Furthermore it can be used as a basis for comparison of the effectiveness of the actual identification search strategy deployed in terms of the accuracy of returned parameter estimates. The problem of inertia (J) and friction (B) parameter extraction of an actual BLMD system over a sinc function shaped multiminima cost surface [19], based on step response feedback current (FC) target data which has a constant amplitude swept frequency characteristic, is investigated as a test case using response surface simulation.

The global minimum estimation, from response surface simulation discussed in section 2.0 below, is targeted towards offline identification of the fixed dynamical load possibly encountered by an embedded BLMD system during the setup and commissioning phase. This is necessary for optimal tuning of the installed BLMD velocity and position loops in any high performance electric vehicle and industrial application. The present work on optimal parameter estimation is mainly concerned with the offline identification of the worst case inertial load that could possibly be experienced by an installed embedded BLMD. This is articulated here through BLMD simulation in torque control mode, using the full reference model developed in [1, 2], and drive experimental step response measurements with three known test cases of shaft load inertia, for validation of the accuracy of the parameter identification strategy, corresponding to:

  1. the no-load rotor inertial value JT =Jm,

  2. medium shaft load inertia JT ~ 4Jm and

  3. large shaft load inertia JT ~ 7Jm.

where JT is the total inertia consisting of rotor Jm and additional shaft load Jl with


The problem of a numerically ‘noisy’ multiminima cost function resulting in non optimal parameter convergence because of local minimum trapping, associated with the adoption of the BLMD reference model in [1, 2] during simulation, in motor parameter identification is examined [20, 21]. An explanation is provided as to the existence of ‘false’ local minima plurality with inaccurate resolution of PWM edge transition times, associated with the choice of fixed step sizes Δt in BLMD model simulation, in both the current feedback Ifj and shaft velocity Vωr MSE objective functions. An explanation is also furnished as to the existence and proliferation of genuine local minima with the observed feedback current (FC) target data used in penalty cost surface generation, which will be shown to posses an inverted sinc function-like shape. Details are presented, through MSE response surface simulation with coarse step sizes chosen initially for the inertia J and friction B parameters employing shaft velocity (parabolic cost surface) and feedback current (sinc-like surface) experimental target data respectively, to shed light on the numerical noise problem for SI purposes. Both simulated MSE response surfaces reveal on a macro-scale the presence of a ‘line minimum’ of possible feasible solutions in a stationary region, enveloping a global extremum within the central surface fold, principally in the B-parameter direction. A novel mathematical approximation [17], which provides verification of the cost surface shape in both cases, is given and is used to provide information on the existence of a unique global minimum with an accompanying optimal parameter set X¯opt={J¯opt,B¯opt}T instead of a multiplicity of candidate options, X¯optj={J¯opt,B¯optj}T, along a ‘B - line minimum’, for j = 1,2 …. Details of BLMD model simulation at a finer parameter step size δX, which illuminates the problem of a noisy cost surface, are also provided for both objective functions. An independent statistical analysis appraisal of the computation ‘noise’ voltage engendered in the search for accurate PWM transitions, based on a novel theoretical estimation [18] for the random error pulse energy expectation associated with PWM replication with chosen simulation step size Δt, is also provided. This probability analysis in itself provides a useful insight into the induced noise mechanism with chosen time step size and highlights the magnitude contribution of the random error ‘noise’ voltage with PWM resolution to the overall accuracy in the BLMD model simulation exercise. The effect of inherent ‘noisy’ evaluation of the PWM edge transition times during BLMD simulation is transferred as a lack of smoothness in the simulated construction of the MSE cost surface at the micro-scale for very low step changes δX in the BLMD dynamical parameters J and B.

A novel mathematical analysis [21] is presented, via embedded quadratic curve fitting in the MSE cost surface, to establish the worst case parameter quantization step size δXL necessary to overcome cost function ‘noisiness’. This analysis also provides a radius of convergence rX in parameter space about the global minimum for any parameter identification search strategy and establishes a bound on the limits of accuracy for the returned optimal parameter estimateX^opt={J^opt,B^opt}T. Furthermore this methodology provides a sensitivity measure of the MSE cost surface selectivity for both the step response shaft velocity and current feedback response surfaces in the neighbourhood of the global extremumX¯opt. This surface variability metric dependency on elemental parameter variation δX can then be used to decide on the best objective function for parameter extraction purposes based on the accuracy of the returned estimate. The choice of the FC target data is explained for its excellent coherence properties, based on frequency and phase attributes from step response tests, in checking BLMD model fidelity and accuracy and also for its high selectivity in penalty cost function formulation for accurate parameter identification. Furthermore it will be shown that there is an improvement in FC cost surface selectivity with longer data training records while the converse effect is manifested for shaft velocity target data with measurement data length in the reduction of cost surface curvature in the vicinity of the global minimum. These current feedback step response attributes arbitrate in its favour as the most suitable choice of target data in MSE cost function formulation.

In the absence of embedded drive application details from the BLMD manufacturer [17] no precise limits on the desired accuracy of the returned J and B parameter estimates could be affixed to the parameter identification strategy for velocity controller tuning purposes in the commissioning phase. However the use of a quantized metric δXL, as mentioned previously, in parameter space puts a limit on the parameter resolution accuracy possible during identification of the BLMD dynamics in electric vehicles. It should be noted that without the imposition of this parameter quantization strategy there is a risk of false minimum trapping of the identification search algorithm in a ‘noisy crevice’ [18] in a side-wall of a cost surface, besides local minimum capture, well away from the global minimum estimate. This novel quantization procedure in parameter space, which eliminates the effect of simulation step size related computation induced noise, results in the availability of a smooth cost surface over which a parameter identification search algorithm will work and converge to an optimal estimate [17,18,21]. One other benefit of the parameter quantization process is that it divides up parameter space and restricts the identification strategy to a countable number of parameter lattice points [18] and thus minimizes the search time to global optimality.

A further aspect of concern besides false minimum trapping is that all optimization algorithms for BLMD parameter identification proceed in a continuous search of parameter space to a convergence estimate of the parameter vector sought with an end stopping criterion [22,23]. The norm of cessation of the optimal parameter search strategy is generally based on the smallness of cost reduction over successive iterations within a specified error bound ε at termination. The termination criteria are generally not focused on the smallest percentage variation of the parameter estimates acceptable. However with the quantization δXL of parameter space for response surface smoothness, limits for parameter resolvability can be imposed by restricting the identification search process to an integral number k of quantum steps kδX commensurate with the percentage accuracy %X required in absolute terms such that %X = kδX. This restricted step approach, in terms of the specified parameter accuracy sought for BLMD tuning purposes during the setup and commissioning phase, can reduce the SI computation time to optimality [18].


2. Response surface simulation and analysis [17]

The concept of a simulated response surface (RS) is presented as an aid to motor dynamical parameter optimization in high performance Brushless Motor Drive (BLMD) identification with a multiminima objective function. This methodology provides useful information concerning the formulation and nature of the most suitable objective function to be minimized, based on actual drive experimental test data available and BLMD model simulation, coupled with an effective system identification (SI) strategy for accurate motor parameter extraction. This simple approach, although computationally intensive, can also be used as an alternative means of providing the optimal set of parameter estimates from inspection of the global minimum location on the simulated cost surface with embedded local minima. Furthermore it can be used as a basis for comparison of the effectiveness of other identification search strategies deployed, such as the Powell Conjugate Direction search method [18] and Fast Simulated Diffusion algorithm [20,21], in terms of the accuracy of returned parameter estimates. The problem of inertia J and viscous friction B parameter extraction of an actual BLMD system over a sinc-function (sinx/x) shaped multiminima cost surface, based on step response feedback current (FC) target data which has a constant amplitude swept frequency characteristic, is investigated using response surface simulation. The choice of the FC target data is based on its excellent coherence properties [24] from step response testing, for checking BLMD model fidelity and accuracy and for the penalty cost function formulation in SI. This difficulty with a multiminima objective function converging to a non optimal parameter estimate, associated with the adoption of the FC target data for motor parameter identification, is examined in the FSD method [20, 21]. An explanation is provided as to the existence of local minima plurality with the observed FC target data used in the sinc-like penalty cost surface generation. All classical optimization techniques [22], with the exception of modern statistical methods [21], are known to have difficulty with this type of cost surface in identifying the optimal parameter vector. The problem arises with initialization of the search strategy far from the global minimum resulting in possible local minimum trapping and non optimal convergence of the cost minimization algorithm during the parameter extraction process. This response surface [RS] methodology, however, provides a simple and effective alternative to classical methods in acquiring an accurate estimate of the global minimum. Results are presented, which demonstrate the efficacy and reliability of the RS method in returning accurate estimates for ‘known’ values of the BLMD shaft dynamics. The application of this FC step response related multiminima cost function in parameter extraction is compared with the alternative parabolic shaped shaft velocity objective function for cost surface selectivity in the vicinity of the global minimum and for accuracy of the returned identified parameter estimates. A mathematical approximation analysis is provided for verification of the cost surface shapes resulting from the deployment of step response FC and shaft velocity as target data in objective function formulation.

2.1. Cost function formulation [18]

Response surface simulation is a useful graphical tool [25] in system identification and can easily be applied to motor parameter extraction and BLMD model validation. This visual concept, which has been used in process control optimization [25], provides an intuitive insight into the topographical structure of the cost function to be minimized and the rapid location of the global minimum. It also provides information on the most suitable identification search strategy that should be adopted in parameter space to obtain an accurate estimate X^opt={x^1 opt,x^2 opt}T of the motor dynamics where the inertia Jx1 and viscous friction Bx2 are the coded variables. The location of the global minimum stationary point can be obtained by inspection from the simulated cost surface. This approach, although computationally expensive, can be used to secure an independent alternative optimal estimate X¯opt={J¯opt,B¯opt}T as a reference against which the accuracy of other parameter identification search schemes such as the Fast Simulated Diffusion [26] can be judged.

BLMD parameter extraction is generally based on the minimization of the errors of fit ek between the observed motor drive target data and BLMD model responses in terms of the controlled parameter vector X. This identification process results in the adjustment of the J and B parameters towards global optimality. The search strategy is performed in the neighbourhood of the global extremum using the least squares error criterion in the cost function formulation between each value of a time series

t={tk|tk=kT, T=Δt, 1kn, kN}E2

of n experimental sample points of the actual motor drive response g(tk) as the target data reference and the corresponding simulated model response f(X, tk). The objective function E(X) is defined as the mean square error (MSE) from the residual vector as






The MSE generates an error response cost surface in parameter X space based on target data from one of the internal test points in [1]. The Powell Conjugate Direction (PCD) [23] and Fast Simulated Diffusion optimization techniques [27] can be applied in conjunction with the BLMD model in [1] to the response surfaces corresponding to motor shaft velocity Vωr and winding FC Ifa target data respectively, obtained in torque control mode for different shaft load inertia listed in [1], for optimal parameter Xopt extraction.

Figure 1.

Experimental FC Cost Surface

Figure 2.

Shaft Velocity Cost Surface

The simulated response surfaces E(X) are derived from BLMD simulation, with a fixed time step of 1µs and appropriate decimation factor, using the model test point o/p f(X,tk) in conjunction with the sampled experimental target data


as the target reference. These penalty cost functions are depicted in Figs.1 and 2 for zero shaft load conditions over parameter space X=[J,B]T with a crude mesh size δX chosen as per Table 1 to initially determine surface shape, according to the rotor inertia and friction tolerances likely to be encountered in practice. The experimental test data training records used in the MSE formulation for each objective function are displayed in Figs.1 and 2.

MSE Cost Surface Type E(X) Current Feedback: Eifa(J,B) Shaft Velocity: Eωr(J,B)
Data Training Record g(tk) Current Feedback: Ifa(tk) Shaft Velocity: Vωr(tk)
No. of Data Points Nd @ 20μs 4095 4095
BLMD Parameter varied x Jm (kg.m2) Bm (Nm/rad/sec) Jm (kg.m2) Bm (Nm/rad/sec)
Nominal Parameter Value xm 2.8x10-4 2.14x10-3 2.8x10-4 2.14x10-3
Parameter Tolerance Band Δ x ±20% ±80% ±20% ±80%
Crude Parameter Step size δ x 1.33% 4% 2% 4%
No. of Parameter Steps Nx 30 40 20 40
Parameter Value Returned 2.99x10-4 ~1.54x10-3 3.024x10-4 ~1.626x10-3
Assumed Optimal Parameter Vector Xo for Response Surface Analysis
Xo 3.0x10-4 2.14x10-3 3.0x10-4 2.14x10-3

Table 1.

Experimental Cost Surface Formulation for Zero Shaft Load (NSL) Conditions

Figure 3.

BLMD Current Feedback

Figure 4.

Rotor Shaft Velocity Vω

The anticipated variation in the search cost, likely to be encountered during BLMD system identification (SI) over the parameter tolerance band of interest, can be gauged from cross sections through the chosen response surface at nominal values of the rotor parameters [Jm,Bm]T. The cost variations associated with specific dynamic parameters are illustrated in Figs.5 and 6 for motor current feedback and in Figs.7 and 8 for shaft velocity target data. These cross sections provide important information regarding the surface shape and curvature and consequently about the nature of the stationary points found and type of SI search algorithm that should be deployed over such hitherto surface ‘terra incognita’.

Figure 5.

MSE-Ifa Cross section at Bm

Figure 6.

MSE-Ifa Cross section at Jm

Figure 7.

MSE-Vωr Cross section at Bm

Figure 8.

MSE-Vωr Cross section at Jm

The FC cost ‘landscape’ highlights the existence of several parabolic shaped ridges, interspersed with embedded synclines within its sinc-like folded topography, with a consequent plurality of local minima. The cost terrain also shows the presence of a stationary elliptical shaped ridge system centrally located in the contour map of Fig. 9 with the possible existence of a ‘line minimum’ [25] along the principal/major axis. These multiminima folds are manifested in the constructive and destructive interference patterns encountered in the frequency ramp up of the FC sinusoid, when compared with the optimal parameter reference or test data waveform, during the transient phase of motor acceleration. The shaft velocity cost surface is parabolic shaped as seen from the contour map in Fig. 10 but is less selective than its FC equivalent in the vicinity of the global minimum when the respective cost surface cross sections with equivalent parameter grid sizes are compared.

Figure 9.

Experimental MSE-Ifa Contour Map

Figure 10.

Experimental MSE-Vωr Contour Map

It is evident from Figs. 9 and 10 that both objective functions possess long wedge shaped stationary valleys in the response surfaces with no ‘apparent’ clearly defined global minimizer. The observed near linear dependence of the surface shape on the parameters in a ‘line minimum’ along the valley floor indicates that B is commensurate with J in the ratio J/B which is the dynamical time constant τm of the motor. The B parameter, which is the least likely of the two to vary in the electromechanical drive applications [8,9] can to be acquired from dynamic testing as per [1] to free the other parameter for identification purposes. This reduces the identification problem to single parameter extraction in J or alternatively in τm, where parameter decoupling is non essential, for controller design purpose.

Response surface simulation provides an alternative route of accurately estimating the optimal parameter vector Xopt by means of inspection of the surface minimum cost. This method, although computationally expensive, can be used as a yardstick by which the overall convergence performance of other identifications schemes [21] can be contrasted, such as FSD, over a range of motor shaft inertial loads. The response surfaces can be simulated initially using a coarse parameter mesh size, for a range of supposedly ‘unknown’ motor inertial load test cases for shaft velocity and current feedback MSE objective functions, for rapid location of the global minimum. Further refinement in mesh size can be made down to the parameter step sizes necessary in the vicinity of the global minimum for accurate resolution of the optimal parameter set. Results, which demonstrate the accuracy and effectiveness of RS simulation, are presented for global minimum estimates of motor shaft inertia which are in close agreement with known test inertial load values.

2.2. Novel mathematical analysis of response surface [18] – Modelling and simulation

Response surfaces can be generated for the BLMD shaft velocity and current feedback step responses, as the mean squared error cost function between an actual drive experimental target data record and simulated model responses, by varying J and B over the two dimensional dynamical parameter space of interest. This graphical procedure is then used to shed light on the shape of the respective cost surfaces and to make a decision as to the most efficient parameter identification strategy to be deployed in each case. Inspection of each of the 2-D MSE response surfaces reveal the existence of ‘open’ wedge shaped stationary regions principally in the B-parameter direction containing what appears to be a global ‘line’ minimum in both cases. From a parameter identification perspective such open stationary regions would mean an infinite number of admissible solutions and thus uncertainty in the parameters extracted. The presence of such a difficulty would require careful measurement of one the parameters, in this case the friction as this is the principal direction that the line minimum appears to exists, in order to free the other (J) for identification. A novel mathematical analysis is presented in this chapter to determine whether or not these embedded stationary regions are open. This approach is articulated by formulating a simple quadratic model approximation of the cost surface stationary regions over a small neighbourhood of parameter space, with interacting J and B terms, for proposed model accuracy. The BLMD model step responses are also approximated by simple analytical expressions over response time spans that are very short by comparison with the dynamical time constant τm for validation and accuracy of the response surface quadratic model approximation. These simple step response representations, in which the parameters J and B can be adjusted over the space of interest for local cost surface generation and analysis of the stationary region, are included along with the relevant experimental target data in the cost surface quadratic model approximation. This mathematical analysis, employing the simplified quadratic model for both cost surfaces, can be used to show:

  • that the stationary regions for the current feedback and shaft velocity objective functions are closed and bounded indicating the presence of a trapped global minimum,

  • how closely the dynamical J and B parameters are coupled by making a comparison of the extracted quadratic model eigenvalues,

  • that a line minimum exists principally in the B parameter direction and quantifies the extent of this B-line minimum by the eigenvalue ratio

  • establishes the degree of ill conditioning for the global minimum solution parameter vector estimate XS extracted from the minimized quadratic model.

Furthermore this analysis also demonstrates that the current feedback response surface has better selectivity in the global stationary region than the shaft velocity equivalent with increasing data record lengths. This outcome helps in the decision analysis that favours the use of current feedback target data in cost function formulation for dynamical parameter identification.

Figure 11.

EM Torque Variation with Bm & Jm

Figure 12.

Simulated Shaft Velocity

The observed topographical features in the above penalty response surfaces can be anticipated from the following approximation analysis. Initially the developed electromagnetic torque Γe is at a maximum for unit torque demand step input Γd and remains so for a very short time as per the BLMD model simulation in Fig.11 until the shaft speed starts to build up exponentially as in Fig. 12 with time constant τm. The back-emf term vej in [1] becomes substantial causing a decrease in winding current ijs which reduces the applied torque. Furthermore the increased rotor angular velocity ωr causes the machine impedance angle ϕz in [1] to approach π/2 and forces the winding currents into quadrature with the current command signals idj with subsequent torque reduction as in [1]. The variation in applied motor torque with the worst case spread of dynamical time constant τm values, observed for the parameter tolerance ranges in Table 1 with zero shaft load conditions, is small over the motor acceleration period (t^=0.08sec60%τm) shown in Fig.11. The average value of applied mechanical torque Γem is 1Nm and is assumed constant over the period t^ for tractability reasons in the following analysis of the cost surfaces used in the PCD and FSD methods of parameter extraction. However this value deteriorates over longer time spans as the winding current moves out of phase alignment with current demand as motor speed increases and thus with the back EMF. The simulated shaft speed variation with time, based on the nominal parameter vector Xm in Table 1 and displayed in Fig.12 for a step i/p torque demand Γd (~1v) is given by (7)-(a)

ωrm(t)=Km(1etτm) with  Km=Γem/Bm and  0tt^ (a)ωr0(t)=K0(1etτ0) with K0=ΓemB0;  τ0=J0B0. (b)E7

Similarly the corresponding shaft speed variation with time at the assumed optimum parameters {Jo, Bo} in Table 1, which are be identified from cost surface trial analysis, is given by (7)-(b)

The sampled motor speed ‘test’ data ωro(tk) generated via (7)-(b) can now be used as target reference ‘test’ data in the simulated trial cost function EωrO for analytical purposes. The optimal parameter set Xo is supposedly unknown and the task here is to obtain a good estimate [J¯o,B¯o]T of this vector in the following cost surface analysis for verification of the RS strategy. The variation in the time constant τm over the permitted parameter tolerance ranges employed in the response surface generation, such as those in Figs 1 and 2 relying on experimental test data, is insufficient to cause departure from nominal applied torque Γem for the short time span shown in Fig.11. The shaft speed variation in this instance is approximated by

ωr(t)=K(1etτ) with K=ΓemB;  τ=JB.E8

Figure 13.

Simulated Velocity Cost Surface

Figure 14.

BLMD Current Feedback

The resulting MSE cost function construct, illustrated in Fig.13 with details in Table 2, is for simulation purposes given by


with target data ωrm. The parabolic cost variations associated with specific dynamic parameters for shaft velocity target data are illustrated in Figs.15 and 16. The corresponding winding current feedback ifa(t) has the characteristics of a frequency modulated sinusoid during the exponential buildup of motor shaft speed in that it exhibits the features of a constant amplitude swept frequency waveform as shown in Fig.14. The effect of shaft speed increase on the phase angle ϕ of the FC response is determined from (8) as


Figure 15.

MSE-Vωr Cross section at Bm

Figure 16.

MSE-Vωr Cross section at Jm

Figure 17.

Simulated FC Cost Surface

The frequency modulated FC, which is current regulated by GI in [1], is given by


with If ≈1 amp for a unit step torque demand i/p. The resultant FC cost surface generated from simulation in Fig.17, with parameter grid sizes in Table 1, is based on the target shaft velocity ωrm(t) in (7)-(a) for nominal values of the dynamical parameters Xm with


Figure 18.

MSE-Ifa Cross section at B*

Figure 19.

MSE-Ifa Cross section at J*

The sinc-profile cost variations associated with specific dynamic parameters for motor current feedback target data at nominal values of the BLMD parameters [Jm,Bm]T are illustrated in Figs.18 and 19. The MSE penalty cost function can described in a more general form about Xm as


with either target data training record deployed using the representation


The nature of the global stationary region embedded in either cost surface, described by (9) or (12), can be explored in canonical form [25] by fitting a quadratic model using a Taylor series. This two dimensional truncated series expansion, with quadratic terms measuring the surface curvature, is anchored at the nominal value Xm to establish the principal axes/directions in parameter space for global minimum search. It is assumed that the expansion pivot Xm is in proximity to the supposed global optimum XO in the case of the FC objective function as this consists of parallel ridges interlaced with folds enveloping local minima regions which obscure the global extremum position. The response surface model ℜ f can be expressed, in either case with (14), in terms of the variables Jx1 & Bx2 and low order interactive terms βij about Xm as

f=β0+β1(x1x1m)+β2(x2x2m)+β112!(x1x1m)2     +β222!(x2x2m)2+β12(x1x1m)(x2x2m)+εE15

with random modelling error ε. The surface model can alternatively be approximated in compact matrix form as


with constant coefficient matrices determined from the cost at Xm, based on target data


by the gradient vector B given by


and the symmetric Hessian matrix G^


which determines the curvature in the vicinity of a local minimum via


The set of constant coefficient differential equations pertaining to (15) are obtained via (13), using either target data record (7)-(b) or (11) with Ifa(t)|ωr(t)=ωr0(t), as


The required first and second order partial differential equations, based on the shaft velocity ωr, to substantiate expressions (22) to (26) are given by

Target Data Record Length Nd with Time Step 20μs 2000 Points - ωr0(tk) 2000 Points - Ifa0(tk)
t^=0.04sec ~31%τm t^=0.04sec ~31%τm
Target Data Parameters
X0 =[J0,B0]T “To be identified
Shaft Velocity Reference Data [3.0x10-4, 2.14x10-3]T Current Feedback Reference Data [3.0x10-4, 2.14x10-3]T
Quadratic Model Fulcrum
Xm =[Jm,Bm]T
Model Surface R^ωr:
[2.8x10-4, 2.14x10-3]T
Model Surface R^Ifa:
[2.8x10-4, 2.14x10-3]T
Model Cost R^mf at Xm 19.553 0.098
Constant β0 via (21) 19.553 0.098
Gradient Vector B [β1, β2]T
via (22/3)
[-2.079x106, -3.279x104]T [-9.827x103, -113.345]T
Hessian Matrix G^ [β11β12β12β22]
via (24/5/6)
[1.238x10111.958x1091.958x1098.873x107] [4.592x1085.114x1065.114x1061.08x105]
Stationary Point
Xs =[Js,Bs]T via (38)
[2.968x10-4, 2.138x10-3]T [3.005x10-4, 2.216x10-3]T
Slope at Xs via (37) [9.313x10-10, 1.455x10-11]T [0, 1.421x10-14]T
Model Cost R^sf at Xs 2.086 -7.654x10-3
Quadratic Form Q(Xs -Xm) via (39) 34.934 0.211
Normal Form of G^
Transformation/Modal Matrix T
with T1G^T=Λ
Normalized Eigenvectors
Normalized Eigenvectors
Co-ordinate Rotation θ -1.813º -1.276º

Table 2.

Summary of Cost Surface Quadratic Modelling Details at Xm

Target Data Record Length Nd
with Time Step 20μs
4095 Points t^=0.082sec ~62.6%τm 4095 Points t^=0.082sec ~56.5%τm
Target Data Parameters
“To be identified”
Shaft Velocity Reference Data
for zero shaft Inertial load (NSL)
Fig. 32; Ref [1] below
Current Feedback Reference Data for zero shaft Inertial load (NSL)
Fig. 29; Ref [1] below
Quadratic Model Fulcrum
Xm =[Jm,Bm]T
Model Surface R^ωr:
[2.8x10-4, 2.14x10-3]T
Model Surface R^Ifa:
[3.1x10-4, 2.14x10-3]T
Model Cost R^mf at Xm 66.543 0.081
Constant β0 via (21) 66.543 0.081
Gradient Vector B [β1, β2]T
via (22/3)
[-6.842x106, -2.422x105]T [1.865x104, 449]T
Hessian Matrix G^ [4.02x10111.376x10101.376x10108.801x108] [1.809x1094.061x1074.061x1072.843x106]
Stationary Point
Xs=[Js,Bs]T via (38)
[2.964x10-4, 2.159x10-3]T [3.00x10-4, 2.124x10-3]T
Slope at Xs via (37) [4.657x10-9,1.746x10-10]T [4.002x10-11, 5.116x10-13]T
Model Cost R^sf at Xs 8.232 -0.015
Quadratic Form
Q(Xs -Xm) via (39)
116.624 0.193
Normal Form of G^
Spectral Condition No. η 0.985x103 0.937x103
Contour sign check (51/2) -1.645x1020 -3.495x1015
Contour Eccentricity e 999.999 x10-3 999.999 x10-3
Modal Matrix T [999.41334.24634.246999.413]103 [999.74822.46122.461999.748]103
Co-ordinate Rotation θ -1.9624º -1.2874º

Table 3.

Details of Cost Surface Quadratic Model Fit at Xm based on actual BLMD Experimental Test Data shown in [1] for Zero Shaft Inertial Load Conditions

The corresponding set of partial derivatives with FC Ifa are obtained via (11) as

2IfaB2=Ifap2(tB[2Kωr]2(τB)ωr)2           sinp(Ktτωr)p(2tB[3K2ωrB]6(τB2)ωrτ(tB)ωrJ)E35
2IfaJB=Ifap2(tB[2Kωr]2(τB)ωr)(ωrBτωrJ)           sinp(Ktτωr)p(2ωrB2+(tB)ωrJ+2(τB)ωrJ)E36

The variation of the directed contour gradient over the fitted cost surface model, given by


is used to locate the global optimum X0 in the parameter hyperspace region of interest. The condition necessary [22] for the presence of a stationary point Xs is the existence of a vanishing gradient in the neighborhood of Xm located within the parameter tolerance band ΔX with


from (37) and the nature of which is determined by the local curvature from the sign of the quadratic form [28]


The parametric details, which include estimates of the gradient vectors and Hessian matrices at Xm for the indicated data record lengths, of the fitted models to the cost surfaces illustrated in Figs.13 and 17 are summarized in Table 2. Similar parametric quantities, employing BLMD experimental test data, are given in Table 3 for cost surface models shown in Figs.1 and 2.

2.2.1. Novel analysis of global minimum estimation and response surface selectivity [18]

An estimate of the cost surface global minimum X^opt is provided in each case by inference from the vanishing gradient in (37) with location of the fitted model stationary point Xs in (38). A sufficient condition for the existence of a global minimizer at Xs is that Q(Xs - Xm) must be positive-definite [28] in which Q(Xs - Xm) > 0 for XsXm. This is verified by the sign of the eigenvalues λi of G^ in Table 2 which are determined from the characteristic equation


The accuracy of global estimates returned in each case for the inertial parameter J in Table 2 admit to the quality and goodness of fit of the models employed for cost surface approximation in the vicinity of the global extremum. The contributory effect of parameter interaction in model approximation in both cases is not insignificant with coefficients βij comparable in magnitude to the geometric mean of the eigenvalues of Hessian G^ in Table 4 defined by


The relative magnitudes of the Hessian curvature components provide information about the uniqueness of the solution Xs in (38), via the matrix condition number in Table 3, based on the infinity norm defined as

cond(G^)=G^G^1 where G^=max1in(j=1n|βij|)E42
Uniqueness of Global Minimum Estimate
Cond G^(Xm) Shaft Velocity Reference Data 2.2126x103 Current Feedback Reference Data 9.1878x103
Spectral Condition No. η 2.144x103 8.99x103
Geometric Mean λ^ 2.674x109 4.844x106
Cost Surface Selectivity and Fitted Model Re-evaluation at Global Minimum Estimate Xs
Fitted Model Fulcrum X s [2.968x10-4, 2.138x10-3]T [3.005x10-4, 2.216x10-3]T
Model Constant β0 at X s 0.376 4.275x10-4
Gradient Vector at X s [3.562x104, 579.78]T [11.535, 0.081]T
Re-evaluation of G^at X s [9.142x10101.437x1091.437x1093.145x107] [4.295x1084.99x1064.99x1066.08x104]
Global Estimate Update Xs1 [2.996x10-4, 2.138x10-3]T [2.998x10-4, 2.159x10-3]T
Slope at Xs1 via (37) [1.717x10-9, 2.547x10-11]T [8.413x10-12, 9.948x10-14]T
Eigenvalues of G^at X s [9.144x1010008.846x106] [4.295x108002.819x103]
Modal Matrix T [999.87615.72315.723999.876]103 [999.93311.61811.618999.933]103
Residual Cost R^sf at Xs1 6.439x10-3 -3.996x10-6
Quadratic Form Q(Xs1 -Xs ) 0.738 8.629x10-4

Table 4.

Results Derived From Cost Surface Quadratic Fit in Table 2

The matrix condition number is much greater than unity in both cases, with the highest value associated with the FC response surface, which indicates a sizeable measure of ill conditioning in the extraction of the global estimate in (38). The curvature component β11 associated with the J-parameter is much greater than that associated with damping B by about three orders of magnitude which indicates greater selectivity of the solution Js along the J axis. This suggests the presence of many potential solutions to (38) along the B-parameter co-ordinate direction due to poorer selectivity or smaller curvature component β22. A more complete interpretation of the nature of the cost surface syncline containing the stationary point region is obtained from the spectral condition number η of G^ [22] as


The relative magnitude η of the eigenvalues indicate that a ‘line minimum’ of potential solutions, which explains the degree of ill conditioning in the global solution estimate, is feasible due to the ‘long’ elliptical shape of the contour map associated with the stationary point zone of convergence in both cases as shown in Figs 9 and 10. The elliptical character of the response surface model in the vicinity of the global minimum estimate can be visualized by a coordinate translation of the parameter axes to Xs as pivot with


resulting in the modified representation from (16) as

R^f=β0+BT(V+XsXm)+12(V+XsXm)TG^(V+XsXm)      =R^sf+12VTG^V (a)where R^sf=β0+BT(XsXm)+12(XsXm)TG^(XsXm) (b)E45

The normalized eigenvectors Ti, associated with the distinct eigenvalues of the symmetric Hessian G^ as


in Table 2, can be used as an orthonormal basis to transform the parameter axes along the principal directions of the elliptical shaped contour system.

Figure 20.

A: Simulated MSE-Ifa Contour Map B: Ifa Contour Map with Canonical Variables

This rotation of co-ordinates, with origin anchored to Xs, is displayed in Figs 20 and 10 for both simulated cost surfaces to eliminate the interactive terms β12 in G^. The X co-ordinate angular displacement θ in Figure 20B can be evaluated from the conical expression [29] as


using (15) with values listed in Tables 2 and 3 which are very small. The rotation can also be deduced from the directional cosines of the unit column vectors [a^z1a^z2] constituting the modal matrix T via


The normal form of the response surface model can be expressed with substitution of the canonical variable

As R^f=R^sf+12ZTΛZ=R^sf+12k=12λkzk2=R^sf+λ12z12+λ22z22E50

with model cost R^sf at the global estimate given in Table 4. The nature of the fitted quadratic model can be deduced from the shape of the embedded contours, in the vicinity of the global minimum estimate Xs, by inference from the sign of the scalar discriminant invariant which pertains to elliptical conic sections [29] as


in the x1x2 frame or

-λ1λ2<0 with λ12= 0E52

in z1z2 normal co-ordinates. These contours are elliptical for both choices of observed BLMD target data with a negative discriminant in Table 3 based on the model cost at nominal parameter value Xm. Consequently the stationary region enveloping the global minimum is encircled and thus bounded by elliptical contours rather than contained within an open wedge shaped response surfaces with no define convergence zone. The degree of elliptical eccentricity e of the trapped stationary zone quantifies the extent of the ‘line minimum’ of global minimum convergence, congruent with the major axis, as notionally illustrated in Fig. 21.

Figure 21.

Elliptical Contour Bounded Stationary Zone

This can be determined from consideration of Fig. 21 by recasting the expression for the normal form of the cost surface model in (50) into that for an elliptical contour, evaluated at the nominal parameter value Xm with origin at Xs, as


for model cost differential


with intrinsic parameters

a2=2ΔR^m,sf/λ1 and b2=2ΔR^m,sf/λ2E55

The eccentricity e, which measures the degree of ‘flatness’ of the oblate model contour R^mf specified at Xm and thus the ‘linear extension’ of the global minimum X0, is given in terms of the lateral displacement c of the elliptical foci from the global estimate Xs relative to the length 2b of the major axis as


The eccentricity of the contours for the model target data in Table 3 is almost unity, as a consequence of the large spectral condition number η in each case, indicating a very flat distended ellipse. The elliptical contours approximates an extended pencil-like global minimum predominantly in the B parameter co-ordinate direction because the inclination angle θ in (47) is less than 2º. This is qualified by the magnitude of the axial ratio (AR) which defines the extent 2b of the ‘line minimum valley’ along the principal direction of the ellipse in relation to its girth 2a, given by the minor axis length, as


This contrast of stationary region ‘feature sizes’ in parameter X-space is readily identified as the spectral condition number η in Table 3 with a ‘line minimum’ extension ratio of about three orders of magnitude for each target data training record used in the MSE cost surface description.

The slope and curvature matrix G^ of the fitted cost model including its associated eigenvalues are re-evaluated at the acquired global estimate Xs as summarized in Table 4 to gauge the response surface selectivity either along the parameter co-ordinate directions or the principal axes of the normal form. The second iterative estimate Xs1, along with the residual costs given in Table 4, is very close to the global minimum target X0 listed in Table 2 for both cost surface models despite the large condition number in each case. A quantifiable measure of the fitted model selectivity in an ill conditioned stationary region, tagged by a large spectral number η, at discerning the global minimum can be obtained from the surface curvature κj along a particular parameter co-ordinate direction xj as



R^fxj<<1  at X=XsE59

or alternatively in normal form as


Figure 22.

Ifa Cost Surface Selectivity

Figure 23.

ωr Cost Surface Selectivity

The degree of model selectivity is three orders of magnitude greater in the case of the motor inertia for actual measured target data employed in both response surface approximations as evidenced from the spectral condition number in Table 3 and in Table 4 for simulated target data trials. Consequently this selectivity margin renders a more accurate estimate in the extracted J-parameter which is mirrored by the arguments leading to the feature size ratio in (57). The cost surface selectivity improves along the principal axes of the normal form when the target data length Nd is extended as indicated by the increased magnitudes of the eigenvalues in Tables 4 and 3. This trend in enhanced J parameter selectivity, which is a measure of the accompanying increase in curvature at the global extremum, is displayed in Figs. 22 and 23 for increasing data record lengths and is a manifestation of the narrowing of the cost surface fold containing the directed ‘line minimum’ principally in the B-parameter direction. The selectivity improvements are greater for increased FC step response data record lengths in Fig. 22 than those for shaft velocity target data in Fig.23. This due to the appearance of more FC cycles with reduced periodicity as motor speed increases demanding a greater degree of fitted model accuracy, with smaller margins of error in terms of frequency and phase coherence at the global minimum value, in the extraction of the optimum parameter vector XO during system identification. The shaft velocity step response by contrast losses its excitation persistence with transient speed decay as it evolves towards steady state conditions with increased data capture time. After a sufficient time elapse the target data transient information, responsible for velocity cost surface folding, is submerged by the steady state onset of maximum motor speed conditions. This irretrievable loss of target velocity signal amplitude variation with time results in a reduction of surface selectivity with parameter variation near the global minimum. These considerations admit to a better choice in the current feedback as a suitable candidate for MSE objective function formulation where accurate parameter extraction is essential during the identification phase of optimal controller design in high performance adaptive BMLD systems for electric vehicle mobility. Furthermore the increasing trend towards motor sensorless control [30] obviates the need for separate rotor position sensors with essential information obtained from the motor signature current via FC sensing at the inverter controller o/p. This adoption of sensorless operation in motor drive systems lends added importance to observed FC data as a suitable target function during parameter identification.


3. Response surface noise and parameter quantization

The computation ‘noise’ inherent in the MSE penalty function construct, based on simulated target data at nominal machine parameter values, is manifested as response surface roughness in parameter space. This is due to model nonlinearities and coarseness of evaluation of the PWM switching instants and results in ‘false’ local minima proliferation in the neighborhood of the global minimizer.

Figure 24.

Noisy Ifa Cost Surface with PRBS I/P

Figure 25.

Noisy Ifa Cost Surface with PRBS I/P

A typical example of this is illustrated in Figs 24 and 25 for BLMD model simulations, with and without inverter turn-on delay δ considered, for small step change variations in the stator winding inductance LS and torque constant Kt parameters. These response surfaces were obtained from BLMD simulation, using FC target data for nominal parameter values as in [1], with a 4095 bit maximal length 2.5 volt bipolar pseudorandom binary sequence (PRBS) input stimulus. The response surface in Fig.24 has a very shallow paraboloidal shape for the small parameter tolerance ranges chosen with a rough noisy texture peppered with local minima in the vicinity of the point-like global minimum. The response surface for simulated FC target data is relatively smooth in the absence of inverter delay turn-on with a point-like singularity at the global minimum as shown in Fig.25. The cost functions pertaining to simulated step response FC Ifa and shaft velocity ωr target data, displayed in Figs.26 and 27 for the dynamic parameters {J,B}, are also noisy with point-like multiminima scattered around the ‘pinhole’ stationary point as in the former case. These surfaces are parabolic for very small tolerance ranges selected near the global minimum as in the main lobe of Fig.1 for the FC corrugated surface.

Figure 26.

Noisy Ifa Cost Surface with Step I/P

Figure 27.

Noisy Velocity Cost Surface with Step I/P

The side elevations of the MSE cost functions in Figs.28 and 29 demonstrate very effectively the fractal landscape with multiminima plurality disposed about the global extremum in the FC case.

Figure 28.

Noisy Ifa Cost Surface Side Elevation

Figure 29.

Noisy Velocity Cost Surface Side Elevation

In the simulated velocity response surface shown in Fig.29 a stationary region exists at zero floor cost with no definite observable global minimum point. An alternative perspective of the minimum stationary regions is provided by the contour maps shown in Figs.30 and 31 for FC and shaft velocity target data respectively. The existence of the point-like global minimum singularity with surface noisiness is clearly evident from the level contours in the FC surface relief map. In the case of the shaft velocity response surface the presence of ‘noisy’ local minima strewn over the ‘river bed’ syncline of the global minimum stationary region is clearly defined by the contour map in Fig.31. The occurrence of ‘noisy’ local minima in the above error surfaces presents a difficulty to any classical optimization method in acquiring the global minimizer where fine parameter resolution is concerned.

A more detailed examination of the effect of inverter delay, achieved through BLMD model simulation without current controller o/p saturation using a 1 volt torque demand step i/p, on the one dimensional MSE response surface in Fig. 32 for very small inductance variation reveal a granulated profile which is less pronounced than that in Fig. 33 with the absence of delay.

Figure 30.

Noisy Local Minima Ifa Contour Map

Figure 31.

Noisy Local Minima ωr Contour Map

Figure 32.

Cost Simulation with Delay δ

Figure 33.

Cost Simulation without Delay

The use of a PWM switch transition time search, based on a single iteration of the regula-falsi method to keep simulation time overhead low, marginally reduces the response error as in Figs.34 and 35. The sensitivity [31] of the error response E with inductance Ls


is very low in all cases and for a ±12% inductance variation gives a change of ΔE=1.5×10-4 in 1.75×10-4 for E. Poor cost surface selectivity will ensue in such cases of low sensitivity over a large parameter tolerance range with possible local minimum convergence if the search process is initiated far from the global minimizer with a noisy cost function.

Figure 34.

BLMD Delay & Reg.-Fal. Method

Figure 35.

Zero Delay & Reg.-Fal. Method

3.1. Novel theoretical estimation of PWM edge transition computation noise [18]

A measure of the computation ‘noise’, induced through inaccurate resolution of the PWM edge transition within a simulation time step Δt, can be ascertained from the associated error in random pulsed energy delivery by the inverter to the stator winding within Δt. Since there is one PWM edge transition every half-switching interval TS/2 of the inverter the expectation in the power delivery error to the stator can be obtained [18], from the error in pulse energy dispatch during the time step interval Δt, as


The expected random voltage vn error associated with inaccurate resolution in PWM inverter switching during BLMD simulation is thus given by

vn=E(P)=UdΔt/TS (a)vn=3101/200=21.92 volts (b)E63

If the chosen simulation time step Δt is 1μs and the inverter switching parameters in [1] are substituted into (63)-(a) the expected uncertainty vn in the inverter output voltage Vjg per phase j can be obtained as (63)-(b)

This value of voltage uncertainty in the inverter output is not insignificant as its magnitude is 7.1% of the inverter HT voltage Ud for a simulation time step size of 1μs. The error can be reduced by decreasing the simulation step size Δt for a more accurate resolution of the pulse edge transition time, once its occurrence has been flagged, or alternatively by means of an accurate search using the regula-falsi method as described in [1].

Figure 36.

Ifa Surface Noise with Inertia Variation

Figure 37.

Ifa Surface Noise Variation with B

The statistical considerations of pulsed energy delivery by the PWM inverter in [18], arising from BLMD simulation with a fixed time step size, illuminates the origin of computation ‘noisiness’ and its subsequent manifestation as cost surface roughness as shown in Figs.36 and 37. For fixed shaft load inertia changes, encountered in motive power applications and electric vehicles, the use of a coarse quantization step size δJ in the inertial parameter variable about the nominal value Jm results in smooth generated and noise-free response surfaces as shown in Figs.1 and 2 for actual FC and shaft velocity target data. However for a sufficiently small step size variation in the inertia Jm and damping Bm a ‘noisy’ cost surface with a proliferation of local minima results in both cases as shown in Figs 36 and 37 for corresponding target test data. The degree of resolution of the parameter step size, that can be obtained and then used in an identification search strategy, depends upon the onset of cost surface irregularity.

3.2. Novel mathematical analysis of quadratic curve fitting to noisy mse cost surface [18,21]

The accuracy of any classical identification scheme used in terms of parameter resolving capabilities can be gauged by fitting a quadratic [26] to the response surfaces in each test case and determining rms deviation of the PWM computation noise related residuals. The quadratic fit employed

Q(xk)=b0+b1xk+b2xk2;         1kNE64

for N steps in the indexed parameter xk as shown in Fig.38 with

xk=xm+(km)δx;      xm{Jm,Bm}E65

is based on an infinitesimal step size δx, in model simulation to reflect response surface roughness, and centred in a tolerance band ±Δxm within indexed range m of the nominal value xm as

m=Δxm/δx E66



The nature of the residual error


associated with the least squares quadratic fit to the various cost functions


for example in Figs.39 and 40 for the FC cost surface, is demonstrated by the autocorrelation (ACR) functions


shown in Figs.41 and 42, which are mainly of the impulse type at zero offset, indicating a white noise-like characteristic.

Figure 38.

Parameter Step Size Resolution

FC Target Data with reference to details in Figs.36 and 37
Parameter varied x Jm Bm
Nominal value xm 3.0375×104kg.m2 2.226×103Nm.rad-1.sec
Parameter Step Size δx 0.01% Jm 0.02% Bm
No. of Steps N 200 200
Tolerance Index m 100 100
b0 96.864 3.286
b1 6.373×105 2.891×103
b2 1.049×109 6.516×105
Residual Error WIfa(x)=EIfa(x)QIfa(x) illustrated in Figs.39 and 40
Standard Deviation σ^ 1.916×104 1.127×104
Mean 3.309×108 3.306×109
Peak Absolute Deviation 5.203×104 3.355×104

Table 5.

FC Cost Function EIfa(x) and Quadratic Fit QIfa(x)

The various cost function details with accompanying quadratic fits and corresponding residuals are summarized in Tables 5 and 6, based on FC and shaft velocity test data respectively, for independent parameter variation in the BLMD shaft inertia and damping factor. The quadratic polynomials fitted to the noisy shaft velocity cost surface sections are displayed in Figs. 43 and 44 with coefficients given in Table 6. The corresponding cost residuals associated with the fitted velocity profiles, which appear to be random, are shown for each of the dynamical parameter variables in Figs. 45 and 46.

Figure 39.

Quadratic Error in J @ Bm

Figure 40.

Quadratic Error in B @ Jm

Figure 41.

Error Autocorrelation in J @ Bm

Figure 42.

Error Autocorrelation in B @ Jm

The white noise-like nature of the error-of-fit in the case of the shaft velocity cost surface sections is demonstrated by the impulse characteristic of the ACR spike functions in Figs.47 and 48. The errors-of-fit can thus be considered as a random entity, with an ACR related noise signature, associated with the BLMD simulation model at very high parameter resolution for each of the observed target data records used in the MSE cost formulation. This manifestation is attributed to some residual uncertainty in the BLMD model simulation of the PWM edge transitions at the comparator o/p with dead time, despite the single iteration cycle of the regula-falsi search, which are magnified in the three phase inverter o/p before being fed to the stator winding.

Shaft Velocity Target Data with reference to details in Figs.43 and 44
Parameter varied x Jm Bm
Nominal value xm 3.0375×104kg.m2 2.14×103Nm.rad-1.sec
Parmeter Step Size δx 0.01% Jm 0.02% Bm
No. of Steps N 370 600
Tolerance Index m 100 100
b0 1.698 0.157
b1 1.079×104 1.098×102
b2 1.747×107 2.465×104
Residual Error Wω(x)=Eω(x)Qω(x) illustrated in Figs 45 and 46
Standard Deviation σ^ 1.028×105 8.974×104
Mean 6.982×1010 1.415×1011
Peak Absolute Deviation 3.476×105 2.871×105

Table 6.

Shaft Velocity Cost Function Eω(x) and Quadratic Fit Qω(x)

Figure 43.

Velocity Cost Noise Variation with J

Figure 44.

Shaft Vel. Surface Noise with B Variation

The effect of lowering the drive model simulation time step Δt, as shown in Figs. 49 and 50 for very small parameter variation in the vicinity of the global singularity, translates into a reduction of the MSE as well as gradual removal of response surface roughness. This tangible decrease in surface roughness with time step size, evident form Fig.51, is measured in terms of the standard deviation of the residual errors associated with various quadratic polynomials fitted to each of the FC cost sections. However the computational effort in terms of CPU time increases in proportion with the decrease in time step size for a given simulation trace length. The requirement for surface noise reduction with the elimination of false local minima plurality has to be balanced with a tradeoff in simulation run time in an attempt to reduce computation costs where BLMD model tractability is an issue in parameter identification and as a simulator in practical applications for performance related prediction of proposed embedded drive systems. A Taylor series expansion of the quadratic fit about the parabolic vertex xopt as


with gradient


can now be used to check the limit of parameter resolution and the “radius” of convergence for worst case conditions [19].

Figure 45.

Quadratic Error in J @ Bm

Figure 46.

Quadratic Error in B @ Jm

Figure 47.

Error Autocorrelation in J @ Bm

Figure 48.

Error Autocorrelation in B @ Jm

At best the smallest parameter threshold step size required in simulation to overcome response surface noise, with rms sample estimate σ^, is determined from that value xks near the global minimum as in Fig.38 such that


Figure 49.

Error of Fit Vs Time Step

Figure 50.

Error of Fit Vs Time Step

Figure 51.

Error Variation with Time Step

3.2.1. Parameter quantization and radius of convergence estimation for system identification

The largest threshold step size estimate can be determined by applying Chebyshev’s theorem for statistical measurements [32] to the response surface noise sample [19, 26]. This theorem indicates that at least the fraction 1(1/h2) of all the residuals Wfk in any sample lie within h standard deviations of the mean μ^ with probability


and for h = 4, which exceeds the tabulated peak absolute deviation in all cases in Tables 5 and 6, is 94%. Thus a measure of the worst case parameter resolution is provided by the inequality


for some large xkL via the quadratic minimiser


in (72) as


The parameter resolution limit in terms of quantization step size δx necessary to overcome cost surface noisiness and local minimum trapping in BLMD parameter identification, which is also a measure of the convergence radius rx about the global minimizer in Fig.38, is given by

Parameter varied x Jm Bm
Minimizer xopt 3.038×104kg.m2 2.219×103Nm.rad-1.sec
Minimizer Offset (mkopt) 0.243 16.199
Threshold Locations kL 72 & 129 24 & 143
Worst Relative Step Size δxLxm 0.281% 1.182%

Table 7.

Quantized Step Sizes for FC Cost Function in Figs. 36 & 37

If measurements are referenced to the nominal value xm at the centre of the parameter tolerance range the relative step sizes


of which there are two pending the sign of the quadratic surd in (77), must be corrected by allowance for the global minimum offset

Parameter varied x Jm Bm
Minimizer xopt 3.007×104kg.m2 2.097×103Nm.rad-1.sec
Minimizer Offset (mkopt) -161.146 -201.448
Threshold Locations kL 210 & 312 212 & 391
Worst Relative Step Size δxLxm 0.505% 1.783%

Table 8.

Quantized Step Sizes for Shaft Velocity Cost Function in Figs.43 & 44

The dynamical parameter threshold step sizes δxL, which are by default the convergence radii measures for reliable global parameter estimation, are tabulated for the response surface cross sections in Tables 7 and 8. The resolution of the motor shaft inertia from the tabulated step sizes, which is the most likely to vary and more essential to identify in high performance applications, is higher when the FC cost function is used instead of the shaft velocity equivalent. Convergence of the inertia parameter estimates to the global minimum is enhanced in the former case with a lower uncertainty due to the smaller step size. The degree of selectivity of the fitted response surfaces with respect to the parameter variability [31] given by


can be determined through the sensitivity coefficient


in the vicinity of the global minimum xopt. This measure can then be usefully employed as a performance index to decide on the best target test data available to use in a motor parameter identification strategy. The sensitivities for 2% parameter variability, greater than the largest threshold step size encountered, of the various fitted surfaces are summarized in Table 9. These sensitivity considerations indicate the suitability of FC test data in the objective function formulation for accuracy in parameter identification.

Parameter varied x Jm Bm
FC Response Surface Sensitivities Figs.36 and 37
Sensitivity 28.88 0.82
Shaft Velocity Response Surface Sensitivities Figs.43 and 44
Sensitivity 0.98 0.07

Table 9.

Response Surface Sensitivity

The above method of parameter quantization, employed to surmount cost surface noise and resultant avoidance of local minimum capture during system identification, reduces the search time in parameter space to global optimality. This is due to the reduction of N-Dimensional parameter space into a finite sized hypercube of countable lattice points NC to be searched, within the imposed parameter tolerance bounds ±Δxm, using an interstitial ‘distance’ equivalent to the step size variability in Tables 7 and 8 as


Figure 52.

Simulated Ifa Cost Surface

Figure 53.

Simulated Shaft Velocity Cost Surface

The application of the tabulated parameter threshold sizes δxL in the objective function simulation results in smooth noise-free response surfaces in the stationary region enclosing the global minimum as displayed in Figs.52 and 53. The degree of accuracy achieved by parameter quantization, with restricted step size during dynamical system identification, in acquiring the global extremum Xopt is determined from the critical values in Tables 7 and 8 as the estimate


The accuracy of the estimate in (84) can be improved with the selection of FC target data because of its greater cost surface sensitivity and smaller relative step size. If the length of the test data record is extended with more data values collected, accompanied by a corresponding transient response decay in the observed variables, the selectivity of the FC response surface improves with genuine local minima proliferation while the parabolic Vωr surface concavity decreases. Thus a more accurate global estimate X¯opt is obtained for reference purposes with increased data record length and improved FC surface selectivity. A suitable identification method can then be applied in conjunction with the BLMD model to the response surfaces corresponding to either motor shaft velocity Vωr or winding FC Ifa target data, obtained in torque mode control for different shaft load inertia, in the parameter search process of the optimal estimate X^opt=[J^opt,B^opt]T. The Powell conjugate direction method [22,23] and FSD [19] parameter extraction techniques can be applied, for example, to the respective Vωr and Ifa cost surfaces to obtain X^opt [18].


4. Conclusions

Response surface simulation has been theoretically investigated and shown to be a useful graphical tool in motor parameter identification with a multiminima objective function and BLMD model validation for electric vehicle systems. This visual concept provides an intuitive insight into the topographical structure of the cost function to be minimized, the location of the global minimum, and the relevant identification search strategy to be adopted in parameter space to obtain an accurate estimate. It also provides an alternative parameter measurement strategy against which the accuracy of other parameter identification search techniques can be judged. A novel mathematical analysis of the competing shaft velocity and current feedback response surfaces, for identification purposes, has revealed the existence of a ‘line’ minimum of possible solutions principally in the B-parameter direction via a comparison of the eigenvalues derived from the quadratic model fit of the global stationary region. This analysis also shows that the global stationary region is closed and bounded by elliptical shaped MSE contours, which guarantees the existence of an optimal parameter vector solution. Furthermore a comparison of the quadratic model eigenvalues, for the competing cost surfaces, illustrates the dominance of the current feedback response selectivity and its acceptance as the most suitable objective function during SI for accurate parameter extraction.

The quantization of parameter space to remove ‘false’ local minima proliferation has been examined and demonstrated to be effective in surmounting cost surface ‘noisiness’ engendered during BLMD simulation, with a finite step size, of the PWM natural sampling process. A probability analysis has shown that the error incurred in resolution of PWM edge transition times during BLMD simulation, which is responsible for cost surface granularity, is dependent on the step size and is manifested as a random error voltage at the PWM inverter output. The effect of cost surface selectivity with choice of target data in MSE penalty cost function formulation, for usage in BLMD parameter identification, has been examined with motor current feedback being the preferred option.



The author wishes to acknowledge

  1. Eolas – The Irish Science and Technology Agency – for research funding.

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


  1. 1. R.A Guinee, Extended Simulation of an Embedded Brushless Motor Drive (BLMD) System for Adjustable Speed Control Inclusive of a Novel Impedance Angle Compensation Technique for Improved Torque Control in Electric Vehicle Propulsion Systems in Electric Vehicles - Modelling and Simulations, ISBN: 978-953-307-477-1, InTech ; 2011.
  2. 2. R.A Guinee, Mathematical Modelling and Simulation of a PWM Inverter Controlled Brushless Motor Drive System from Physical Principles for Electric Vehicle Propulsion Applications in Electric Vehicles - Modelling and Simulations, ISBN: 978-953-307-477-1, InTech ; 2011.
  3. 3. Miller, J.; (2010). Propulsion Systems for Hybrid Vehicles, IET, Renewable Energy, 2nd Edition.
  4. 4. R.M. Crowder, Electric Drives and their Controls, 1995, Clarendon Press, Oxford.
  5. 5. A.J. Critchlow, Introduction to Robotics, Macmillan Pub. Co. NY, 1985.
  6. 6. H. Gross, Electrical Feed Drives for Machine Tools, 1983 by Siemens, J. Wiley & Sons.
  7. 7. Moog Brushless Technology User Manual:D31X-XXX Motors,T158-01X Controllers,T157-001 Power Supply, Moog GmbH, D-7030 Böblingen, Feb 1989.
  8. 8. H. Asada and K. Youcef-Toumi, Direct-Drive Robots Theory and Practice, 1987, MIT Press.
  9. 9. R.P. Paul, Robot Manipulators: Mathematics, Programming and Control, The MIT Press, Camb, Mass, USA, 1986.
  10. 10. W.E. Snyder, Industrial Robots: Computer Interfacing and Control, PHI, 1985
  11. 11. M.A. El-Sharkawi and S. Weerasooriya, “Development and Implementation on Self-Tuning Tracking Controllers for DC Motors”, IEEE Trans. on Energy Conv., Vol. 5, No. 1, Mar 1990.
  12. 12. N.A. Demerdash, T.W. Nehl and E. Maslowski, "Dynamic modelling of brushless dc motors in electric propulsion and electromechanical actuation by digital techniques", IEEE/IAS Conf. Rec. CH1575-0/80/0000-0570, pp. 570-579, 1980.
  13. 13. H. Dohmeki and M. Nasu, “ Development of a Brushless DC Motor for Incremental Motion Systems”, Proc. 14th IMCSD annual Symp., pp.63-71, 1985
  14. 14. J.Y.S. Luh, “Conventional Controller Design for Industrial Robots – A Tutorial”, IEEE Trans. On Systems, Man, and Cybernetics, Vol. SMC-13, No. 3, May/June 1983.
  15. 15. K.J. Astrom and T. Hagglund, Automatic Tuning of PID Controllers, Instr. Soc. Amer, 1988, ISBN 1-55617-081-5.
  16. 16. Moog Brushless Technology:Brushless Servodrives User Manual D310.01.03 En/De/It 01.88, Moog GmbH, D-7030 Böblingen, Germany, 1988.
  17. 17. R.A. Guinee and C. Lyden, “Motor Parameter Identification using Response Surface Simulation and Analysis”, Proc. of American Control Conference, ACC-2001,June 25-27, 2001, Arlington, VA, USA.
  18. 18. Guinee, R.A., “Response Surface Methodology”, Modelling, Simulation, and Parameter Identification of a Permanent Magnet Brushless Motor Drive System, Chapter 3, pages 125 – 206, Ph. D. Thesis, 2003, National University of Ireland – University College Cork.
  19. 19. R.A. Guinee and C. Lyden, “A Novel Application of the Fast Simulated Diffusion Algorithm for Dynamical Parameter Identification of Brushless Motor Drive Systems”. IEEE-ISCAS 2000, The 2000 IEEE International Symposium on Circuits and Systems, May 28-31, Geneva, Switzerland.
  20. 20. R.A. Guinee and C. Lyden, “A Novel Application of the Fast Simulated Diffusion Optimization Technique for Brushless Motor Parameter Extraction” UKACC International Conference on Control, Cambridge Univ., Sep 2000.
  21. 21. R.A. Guinee and C. Lyden, “Parameter Identification of a Brushless Motor Drive System using a Modified Version of the Fast Simulated Diffusion Algorithm”, Proc. of American Control Conference – IEEE ACC-1999, San Diego, June 1999, pp.3467-3471
  22. 22. R. Fletcher, Practical Methods of Optimization, 2nd edition,1993, J.Wiley & Sons.
  23. 23. W.H. Press, B.F. Flannery, S.A. Teukolsky and W.T. Vetterling, Numerical Recipes in C, 1990, CUP.
  24. 24. Guinee and C. Lyden, “Accurate Modelling And Simulation Of A DC Brushless Motor Drive System For High Performance Industrial Applications”, IEEE ISCAS’99 - IEEE International Symposium on Circuits and Systems, May/June 1999, Orlando, Florida
  25. 25. R.H. Myers and D.C. Montgomery, Response Surface Methodology - Process and Product Optimization Using Designed Experiments, 1995, J.Wiley & Sons, NY.
  26. 26. R.A. Guinee and C. Lyden, “A Novel Application of the Fast Simulated Diffusion Algorithm in Brushless Motor Parameter Identification”, The 3rd IEEE European Workshop on Computer-Intensive Methods in Control and Data Processing, Sep 7-9, Prague, Czech Republic.
  27. 27. R.A. Guinee and C. Lyden, “Parameter Identification of a Motor Drive using a Modified Fast Simulated Diffusion Algorithm”, Proc. of the IASTED Intern. Conf. on Modelling and Simulation, pp 224-228, May. 1998, Pittsburgh, Pa., USA.
  28. 28. H.A. Taha, Operations Research, Macmillan Publishing Co., NY., 1971.
  29. 29. Protter and Murray, Calculus and Analytic Geometry
  30. 30. J. Holtz, “Sensorless Position Control of Induction Motors - an Emerging Technology, invited paper, IEEE-IECON’98, Proc. of the 24th Annual Conf. of the IEEE Indus. Electronics Society, Aug 31 - Sep 4, 1998, Aachen, Germany.
  31. 31. G. Daryanani, Principles of Active Networks Synthesis and Design, 1976, J. Wiley & Sons.
  32. 32. F. Mosteller, R.E.K. Rourke and G.B. Thomas, Probability with Statistical Applications, 1961, Addison-Wesley Publ. Co.

Written By

Richard A. Guinee

Submitted: 23 August 2012 Published: 19 December 2012