Mathematics and Computers in Simulation 71 (2006) 446–459

Modelling and identification of a non-linear saturated magnetic

circuit: Theoretical study and experimental results
Sandrine Moreau ∗ , Jean-Claude Trigeassou
Laboratoire d’Automatique et d’Informatique Industrielle, Ecole Supérieure d’Ingénieurs de Poitiers, 40,
avenue du Recteur Pineau, 86022 Poitiers Cedex, France
Available online 17 April 2006

Abstract
This paper deals with the phenomenon of magnetic saturation in electrical machines. Taking magnetic saturation into account
leads to non-linear models which need to develop parameter estimation techniques well adapted. At first, the study concerns an
elementary model of saturated iron core coils and the physical parameter estimation of this non-linear model, thanks to a Non-linear
Programming algorithm.
© 2006 IMACS. Published by Elsevier B.V. All rights reserved.

Keywords: Magnetic saturation; Identification; Non-linear Programming algorithm; Magnetic circuit

1. Introduction

In conventional models used in electrical engineering to represent electrical machines, it is quite a common thing
to neglect the magnetic saturation phenomenon, which allows to make the simplifying assumption that inductances
depend on the value of the flux imposed by currents, thanks to the following non-linear relation φ(i) = L(i)·i, where
φ(i) is a picture of the non-linear magnetization characteristic B(H) (Fig. 1) of the employed ferromagnetic material.
Moreover taking into account the saturation phenomenon leads to differential equations which coefficients are time-
varying and consequently to fundamentally non-linear in respect to parameters models. The main difficulty is also to
find a non-linear model structure suited to the studied system and of a complexity level straightforward enough so as
to develop an identification procedure.
As for the parameter estimation, a Non-linear Programming algorithm is imperative.
Firstly, before extending this study to swivel alternating current machines, we contemplate analysing the saturation
phenomenon in the primary winding of a single-phase transformer.
The last mentioned is also represented firstly by a non-linear inductance connected in series with a resistance.
However, it will be seen later that the identification results obtained in the experimental case can be improved if
iron losses are not neglected and taken into account by a resistance in parallel with the inductance. First of all, the
emphasis is put on the mathematical model chosen to represent the magnetic saturation of the inductance.
After developing an iterative identification procedure, both simulation and experimental results are given to validate
the minimization algorithm as well as the model of the saturated inductance.

∗ Corresponding author.
E-mail address: [email protected] (S. Moreau).

0378-4754/$32.00 © 2006 IMACS. Published by Elsevier B.V. All rights reserved.

doi:10.1016/j.matcom.2006.02.006
 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459 447

Fig. 1. Non-linear magnetization characteristic B(H).

2. Saturated inductances modelling

In this study, the magnetizing flux will be described by a simple magnetization curve like Fig. 1 and not a hysteresis
cycle [6]. Thus, iron losses through hysteresis as well as eddy current will be at first neglected. On the same way the
leakage flux is neglected.
A grey-box approach has been carried out to determine a set of analytical functions which approximate at best the
magnetization characteristic B(H) (or the curve φ(i)).
The example of an iron core coil fed by a sinusoidal voltage u(t) and crossed by a current i(t) is considered afterwards.
Fig. 2 represents the primary winding of a transformer with N turns of resistance, R where l and S are, respectively,
the mean length and the section of the magnetic circuit.
Fig. 3 represents the electrical equivalent scheme of the iron core coil described above.
It can be shown from the relation given below that the larger the air-gap e is, the lower is the value of the inductance:
N2
L= (1)
(1/µfer ) + (2 · e/µ0 )
with µfer = µ0 ·µr , where µr is the relative permeability which depends on the induction B.
µ0 is the vacuum permeability (µ0 = 4π10−7 H/m).
With the previous simplifying assumptions, the system is described by the following differential equation:
dφ(i)
U(t) = R · i(t) + (2)
dt

Fig. 2. Scheme of an iron core coil.

Fig. 3. Electrical equivalent scheme of the iron core coil.

448 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459

Fig. 4. Magnetizing curve φ(i).

This equation can be written again:

dφ(i) di
U(t) = R · i(t) + · (3)
dt dt
The term dφ(i)/di = Ldyn (i) represents the dynamical inductance. So Ldyn (i) represents the first derivative of the mag-
netization curve φ(i) (Fig. 4)as it can be seen in Fig. 5.
After several attempts and comparisons with other models [1,2,4,3,7], the mathematical model of the dynamical
inductance has been chosen as follows:
1
Ldyn (i) = c + (4)
a + b · exp(α · i2 )
where a, b, c and α are constant parameters to determine and i is the current in the coil.

Fig. 5. Dynamical inductance in function of the magnetizing current.

 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459 449

i
Ldyn (i) confirms of course the relation φ(i) = 0 Ldyn (µ) dµ.
Ldyn (i) is introduced in Eq. (2) and the following differential equation is obtained:
di(t) R 1
=− · i(t) + · U(t) (5)
dt Ldyn (i) Ldyn (i)
The differential Eq. (5) depends therefore on parameters θ 1 = R, θ 2 = a, θ 3 = b, θ 4 = c and θ 5 = α and consequently
the output model of the current ı̂(t), obtained by solving (5), is fundamentally non-linear in respect to parameters.

3. Parameters estimation algorithm

As the model is non-linear in respect to its parameters, the technique employed to estimate the vector θ- relies on
the minimization of a quadratic criterion thanks to a Non-linear Programming algorithm.
This Non-linear Programming algorithm is the Marquardt one [5] and realizes automatically the compromise
between a Gradient type algorithm far from the optimum so as to get closer to it without instability and a Newton type
algorithm in order to increase the convergence speed in the optimum neighbourhood.
So the minimization of the following quadratic criterion is realized:
K
 K

J= ε2k = (yk∗ − ŷk )2 (6)
k=1 k=1

where εk represents the output error, yk∗ the output measurement and ŷk represents the estimated output.
At the iteration (i + 1), the Marquardt algorithm can be expressed as:
θ- i+1 = θ- i + θ- (7)
where θ- is solution of
 

{ µ · I + Jθθ · θ- = −Jθ }θ (8)
i

with
K

(i∗ (k) − î(θ- i , k))
2
J(θ- i ) =
k=1
K
(9)
∂J(θ- i ) 
Jθ = = −2 · εk · σ- k
∂θ-
k=1

knowing that σ- k = ∂îk/∂θ- is the sensitivity functions vector.

Then the pseudo Hessien of the criterion J is computed in the approximation of Gauss–Newton.
K
 ∂2 J(θ- i ) ∼
Jθθ = = 2 · σ- k · σ- Tk (10)
∂θ- 2 k=1

if µ → 0, a Newton type behaviour is obtained; if µ  1, a Gradient type behaviour is obtained.

The search direction set by µ is determined at each step so as to assure the convergence.
Although the study is here limited to the estimation of only two physical parameters: the resistance R and the
dynamical inductance Ldyn (i), we have to estimate five parameters, which requires to compute five sensitivity functions.
⎤ ⎡
θ1 = R
⎢ ⎥
⎢ θ2 = a ⎥
⎢ ⎥
These five parameters are θ- = ⎢ ⎥
⎢ θ3 = b ⎥ .
⎢ ⎥
⎣ θ4 = c ⎦
θ5 = α
450 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459

The five sensitivity functions are deduced from the general differential equation of the current.
For instance, let us consider the derivative of the differential Eq. (5) in respect to the parameter θ 1 .The sensitivity
function calculated in respect to the parameter R can be obtained by solving the new differential equation given
below:
   
R d(Ldyn (i, θ- )−1 ) ∂(Ldyn (i, θ- )−1 )
σ̇R = − + (u − Ri ) σR + (u − Ri ) (11)
Ldyn (i, θ- ) di ∂R
 
=0

The other sensitivity functions σ- θ = ∂î/∂θ- calculated in respect to the other parameters a, b, c and α, are computed
-
in the same way.

4. Simulation results

The digital simulation of the system constituted by iron core coils fed by a sinusoidal voltage 250 V/50 Hz is carried
out by a fourth-order Runge–Kutta algorithm.
Fig. 6 shows the shape of the simulated current i(t) for a working in saturation mode.
As for the identification protocol, the excitation of the system is achieved by an additional random disturbance, such
as a Pseudo Random Binary Sequence (P.R.B.S.), applied to the amplitude of the feeding voltage.
Moreover, the main property of this P.R.B.S. is to be random in amplitude, this amplitude changing randomly on
the following limited interval: (−50 V, +50 V) in order to excite the nonlinearity of the system.
The parameters of the coil have been set to the following values:
R = 2Ω
1 (12)
Ldyn (i) = 0.03 + H
3.44 + 0.1 · exp(0.3 · i2 )
The sampling period is Te = 0.2 ms.
Identification results are given in the following tables.
Table 1 gives the results obtained by the Marquardt algorithm.
Figs. 7 and 8 represent the estimation convergence of each parameter.
A comparison is made on Figs. 9 and 10 between the curves L̂dyn (i) and φ̂(i) drawn from the estimated parameters
and the accurate ones.

Fig. 6. Shape of the simulated current i(t) in saturation regime.

 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459 451

Table 1
Parameters estimation results obtained in simulation by Marquardt algorithm
Parameters Initial values Estimated values

R (
) 4.1 2
a 3.784 3.44
b 0.15 0.1
c 0.0375 0.03
α 0.45 0.3

Fig. 7. Parameters estimation results.

Fig. 8. Parameters estimation results.

452 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459

Fig. 9. Real and estimated dynamical inductance in function of the magnetizing current.

Fig. 10. Real and estimated magnetizing flux in function of the magnetizing current.

As the estimated curves are very close to the real ones, we can conclude that this simulation study presents a solution
for modelling and estimating non-linear inductances which appear in the saturation of ferromagnetic materials.
And an identification procedure has been proposed and validated in simulation for the parameter estimation of this
nonlinearity linked to the saturation phenomenon.

5. Experimental results

This section puts the emphasis on some practical aspects with regard to the identification of a dynamical continuous
model and presents the experimental set-up which has been carried out on an iron core coil so as to validate the
dynamical inductance model given in Section 2 and the identification algorithm described previously in Section 3.
The experiment is made up of a measurement and formatting system as well as a signals reading–recording system
(Fig. 11).
The first one is composed of current and voltage sensors.
All these measurements are sent to the P.C. via an interface card (Fastlab) which reads and records
the data. The latter contains an analog/digital converter with a resolution of 12 bits and with a conversion
 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459 453

Fig. 11. Experimental system.

time of 1.25 ␮s (800 kHz), a double input–output interface PIO 82C55 associated with a timer 82C52 and a
microcomputer.
Finally signals are converted on 12 bits and transferred to the P.C. via the input–output bus of 16 bits before being
saved and stored on the hard-disk of the P.C. so as to be treated by the identification procedure.
First of all, so as to know approximately parameters values, their values are deduced from power measurements so
as to make at the end a comparison with the results obtained with the identification procedure.
L is determined from the direct knowledge of the reactive power Q, the supplying voltage U and the frequency f
(f = 50 Hz) and R is measured with the volt–amperometric method.
U2
L= (13)
2π · f · Q
We can notice that for the linear part of the magnetization characteristic Φ = f(i), the value of the dynamical inductance
Ldyn (i) = dφ(i)/di is equal to the inductance L(i) = φ(i)/i, since the beginning of the curve can be considered as a straight
line.
To excite the nonlinearity, the amplitude of the feeding voltage is made vary rapidly and randomly on a large
amplitude interval.
For information, the resistance and inductance values of the coil without the iron core are given above (Table 2).
Then power measurements results are given for different thicknesses e1 and e2 of the air-gap with e1 < e2 of the iron
core coil. Tables 3 and 4 sum up the inductance values obtained through power measurements with different air-gap
and voltage values.
Then the results obtained with the parameters estimation method developed above are presented in Table 5.
From the estimation results given above, we can notice that the value of the estimated resistance R is higher than
the resistance value obtained through experimental measurements with the volt–amperometric method.
These estimation results show that it remains an important modelling error because the model based on a saturated
inductance in series with the resistance R does not take into account all the physical phenomena which occur in the
iron core coil. Indeed,

Table 2
Parameter values of the experimental bench
Parameter Data Unit

Resistance (R) 2.7

Inductance (L) 0.013 H

454 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459

Table 3
Power measurements for an air-gap of e = e1
Voltage (V) Reactive power (VAR) Inductance (H)

120 61 0.751
140.5 83 0.757
160 109 0.747

Table 4
Power measurements for an air-gap of e = e2
Voltage (V) Reactive power (VAR) Inductance (H)

120 131 0.349

140 175 0.356
160 229 0.356

Table 5
Identification results for an air-gap e = e1
Parameters Estimation results

R (
) 35.42
a 0.3372
b 1.6129
c 0.3085
α 1.4243

- Ldyn (i) = f (θ- , i) does not represent perfectly the saturation phenomenon,
- the hysteresis is not modelled,
- and particularly iron losses due to eddy currents are not modelled.

So the imperfect modelling leads to an increase of the value of the resistance R, which tries to explain especially
the iron losses firstly neglected in the model.
A solution to this problem consists in improving the model of the iron core coil but a compromise has to be found
so as to keep the model complexity reasonable in respect to the identification procedure.
So as to take into account the iron losses, a second model is developed and used for the identification. This model
presented in Fig. 12 contains a resistance R0 which is added in parallel with the magnetizing inductance.
The identification results obtained with this new model are given below in Table 6.
It turns out that results are improved if the iron losses resistance R0 in parallel with the dynamical inductance is
taken into account.
Now, we can notice that the estimated value of the resistance R is smaller than the estimated one by using the first
model which does not take into account the iron losses.
As the estimated current fits very well with the measured one in Fig. 13, we can conclude that the model of the
dynamical inductance is rather well suited.
Figs. 14 and 15 represent the estimation convergence of each parameter

Fig. 12. Electrical equivalent scheme of the iron core coil with a resistance R0 in parallel with Ldyn .
 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459 455

Table 6
Identification results for an air-gap e = e1
Parameters Estimation results

R (
) 6.46
R0 (
) 2000.09
a −1.015
b 2.38
c 0.051
α 0.33

Fig. 13. Shape of the measured and estimated current i(t) in saturation regime.

Fig. 14. Parameters estimation results.

456 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459

Fig. 15. Parameters estimation results.

As for the parameters a, b, c and α, it is difficult to know in practice if they are well estimated or not because they
have no physical meaning. What matters more is the finale shape of the dynamical inductance obtained with these four
parameters in function of the magnetizing current.
Fig. 16 shows that the estimated value of the dynamical inductance obtained for low magnetizing current – Ldyn is
about 0.76 H – is very close to the value obtained by power measurements (L = 0.75 H).
For the air-gap e = e2 , we can notice from Fig. 18 that the magnetizing current is less saturated than in the previous
case where e = e1 (Fig. 13).
For e = e2 , the estimation of the resistance R becomes equal to its right value, which is of 2.7
(Table 7). As the
magnetic circuit of the iron core coil is less saturated, the modelling error seems to be smaller.
Once more, the estimation results show that the current is rather well estimated, as well as the dynamical inductance
(Fig. 19). Indeed its initial value at low magnetizing current – Ldyn is about 0.36 H – fits rather well with the value of
the inductance obtained through power measurements in linear case without saturation (L = 0.35 H). Figs. 17 and 20
represent the estimated magnetizing flux in function of the current respectively for e = e1 and e = e2 .

Fig. 16. Estimated dynamical inductance in function of the magnetizing current.

 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459 457

Fig. 17. Estimated magnetizing flux in function of the magnetizing current.

Fig. 18. Shape of the measured and estimated current i(t) in saturation regime.

Table 7
Identification results for an air-gap e = e2
Parameters Estimation results

R (
) 2.709
R0 (
) 2000
a −0.558
b 3.043
c −0.01
α 0.06
458 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459

Fig. 19. Estimated dynamical inductance in function of the magnetizing current.

Fig. 20. Estimated magnetizing flux in function of the magnetizing current.

6. Conclusion

In this paper, a solution to model and identify dynamical inductances which vary widely with the magnetic saturation
level is proposed.
A sensible choice of the dynamical inductance model along with a well-adapted estimation technique allow to know
with a better accuracy the behaviour of winding inductances.
This study will be extended to more complicated electrical systems such as the induction motor even if a new
difficulty appears in this case since the number of parameters to estimate becomes greater.

References

[1] J. Faiz, A.R. Seifi, Dynamic analysis of induction motors with saturable inductances, Electr. Power Syst. Res. 34 (3) (1995) 205–210.
[2] A. Keyhani, H. Tsai, IGSPICE simulation of induction machines with saturable inductances, IEEE Trans. Ind. Appl. 4 (1) (1989) 180–189.
[3] E. Levi, A Novel Saturated Induction Machine Model with Application in Analysis of Capacitor Braking, Saint-Nazaire, Electrimacs’96,
September 17–19, 1996, pp. 977–982.
 S. Moreau, J.-C. Trigeassou / Mathematics and Computers in Simulation 71 (2006) 446–459 459

[4] E. Levi, V. Vuckovic, Magnetizing curve representation methods in digital simulation of induction machine dynamics, Int. J. Model. Simul. 10
(2) (1990) 52–56.
[5] D.W. Marquardt, An algorithm for least squares estimation of nonlinear parameters, J. Soc. Ind. Appl. Math. 11 (2) (1963) 431–441.
[6] S. Moreau, Contribution à la modélisation et à l’estimation paramétrique des machines électriques à courant alternatif: application au diagnostic’,
Ph.D., University of Poitiers, 1999.
[7] V. Subrahmanyam, Representation of magnetisation curves by exponential series, Proc. IEE 121 (July) (1974) 707–708.