IN
DEGREE PROJECT ELECTRICAL ENGINEERING, SECOND CYCLE, 30 CREDITS
STOCKHOLM SWEDEN 2017,
Real-time implementation of PMSM software model on external
hardware
ALEXANDER SJÖBERG
Real-time implementation of PMSM software model on external hardware
ALEXANDER SJ ¨OBERG
Master of Science Thesis in Electrical Energy Conversion at the School of Electrical Engineering
KTH Royal Institute of Technology Stockholm, Sweden, May 2017.
Supervisor at KTH: Oskar Wallmark
Supervisors at Inmotion Technologies: Jesper Adolfsson, Per Johansson Examiner: Oskar Wallmark
TRITA-EE 2017:106
Real-time implementation of PMSM software model on external hardware ALEXANDER SJ ¨OBERG
ALEXANDER SJ ¨OBERG, 2017.c
School of Electrical Engineering
Department of Electric Power and Energy Systems Kungliga Tekniska h¨ogskolan
SE–100 44 Stockholm Sweden
Abstract
When developing three phase motor drives, the best way to validate the desired functionality is to connect the inverter to an actual electrical motor. However, when de- veloping functions which are not directly involved in controlling the motor, it could be more efficient to use a real-time software model of the motor. In this master thesis, the de- velopment and implementation of a software model of a permanent magnet synchronous motor (PMSM) is presented. This model was based on general dynamic equations for PMSM in a rotating reference frame (dq-frame). The model was simulated and converted to C code using model based software development in Mathworks Simulink. To provide more realistic performance of the model, a finite element analysis (FEA) was done of an actual PMSM using the software tool FEMM. This analysis resulted in data describing the relation between flux linkage and current which, when added into to software model, limits the produced torque due to magnetic saturation. Both the FEMM model and the final software model was compared to a corresponding actual motor for validation and performance testing. All this resulted in a fully functional software model which was ex- ecutable on the inverter. In the comparison of FEMM model to the real motor, a deviation in produced torque was discovered. This led to the conclusion that the model needed to be improved to perform more alike the real motor. However, for this application the model was considered good enough to be used in future software development projects.
Key words: Model based software development, permanent magnet synchronous motor modelling, real-time software model.
Sammanfattning
N¨ar kontrollsystem till trefasmotorer utvecklas s˚a ¨ar det mest vanliga och troligen det b¨asta s¨attet f¨or funktionsvalidering att k¨ora drivenheten kopplad mot en riktig elektrisk motor. D¨aremot, om funktioner som ej ¨ar direkt kopplade till sj¨alva drivningen av motorn utvecklas, s˚a kan det vara mer effektivt att ist¨allet anv¨anda en mjukvarumodell. I det h¨ar examensarbetet s˚a presenteras en mjukvarumodell av en permanentmagnetiserad synkron- motor (PMSM). Modellen baserades p˚a de generella ekvationerna f¨or PMSM och simuler- ades samt kodgenererades i Mathworks verktyg Simulink. F¨or att g¨ora modellen mer real- istisk s˚a kompletterades den med data som beskriver relationen mellan det l¨ankade fl¨odet och str¨om f¨or att ¨aven ta h¨ansyn till magnetisk m¨attnad. Den informationen simulerades fram i verktyget FEMM genom fl¨odesber¨akningar p˚a en specifik motor typ. Samma motor typ har ocks˚a j¨amf¨orts med den slutgiltiga mjukvarumodellen med avseende p˚a utveck- lat vridmoment vilket resulterade i n˚agot st¨orre skillnader ¨an f¨orv¨antat. Slutsatsen blev s˚aledes att modellen beh¨over f¨orb¨attras f¨or att p˚a ett b¨attre s¨att st¨amma ¨overens med verk- ligheten men att den fungerar tillr¨ackligt bra f¨or den ¨amnade applikationen.
Nyckelord: Modellbasserad mjukvaruutveckling, modellering av permanentmagnetiserade synkronmotorer, realtidsmodellering.
Acknowledgements
I would like to thank my supervisor and examiner, Professor Oskar Wallmark, for all of his support throughout this project. I would also like to thank the staff at Inmotion Technologies and especially Per Johansson, Shuang Zhao, Jesper Adolfsson , Christian Sandberg and Lars L˚aftman for all their involvement in this project.
Alexander Sj¨oberg Stockholm, Sweden August 2017
Contents
Abstract i
Sammanfattning iii
Acknowledgements v
Contents vii
1 Introduction 1
1.1 Background . . . 1
1.2 Project Goals . . . 2
1.3 System set-up and requirements . . . 2
1.3.1 The control system . . . 2
1.3.2 Hardware requirements . . . 3
1.4 Thesis outline . . . 3
2 Theory 5 2.1 Motor equations . . . 5
2.2 Complement with magnetic saturation . . . 6
2.3 Adding flux linkage harmonics . . . 7
2.4 Mechanical model . . . 7
3 Finite Element Method Analysis 9 3.1 The geometry . . . 9
3.2 Extrapolation of BH-curve . . . 10
3.3 Validation . . . 11
3.4 Data extraction . . . 12
4 Simulink 15 4.1 The continuous model . . . 15
4.2 Discretization methods . . . 20
4.3 Fixed point conversion . . . 21
Contents
4.4 Code generation and validation . . . 21
5 Results 25 6 Discussion and Conclusions 27 6.1 Discussion . . . 27
6.1.1 FEMM . . . 27
6.1.2 Simulink . . . 28
6.2 Conclusions . . . 28
6.3 Suggestions on future work . . . 28
6.3.1 Improvement of the FEMM-model . . . 28
6.3.2 Develop the mechanical model . . . 28
6.3.3 Modelling of other motor types (IM) . . . 29
6.3.4 Study the fixed-point conversion work flow . . . 29
A Simulink models 31 A.0.1 Continuous Model . . . 32
A.0.2 Discrete Model . . . 34
B Generated C code 35 B.0.1 Main file . . . 35
B.0.2 Flux linkage to current tables . . . 50
B.0.3 Global defines . . . 52
References 57
Chapter 1 Introduction
1.1 Background
Inmotion Technologies AB designs and manufactures three-phase inverters for electric vehicles like commercial hybrids and forklift trucks. During the development and testing of the inverter, it is necessary to interface and run actual electrical motors in order to vali- date the specific application need and motor parametrization. Interfacing electrical motors requires a mechanical test bench, and in some cases, high currents and voltages. Also, it is not always possible to simulate the, e.g., final vehicle environment in a test bench. In order to simplify and speed-up the design and testing, it would be beneficial to have a parametrized simulation model in the inverter of an electrical motor and its mechanical surroundings. Other benefits with a parametrized simulation model includes that electrical motor parameters can be easily changed for different test cases and that changing vehicle environments also can be simulated. The motor types of interest is first of all induction machines (IMS). IMS is a very popular machine to use in the automotive industry due to its fairly simple design making it easy to build, which also reflects the price. The other motor type is the permanent-magnet synchronous machine (PMSM) which, when com- pared to IMS, are more expensive to build. The reason for this is the permanent magnets which are expensive and also makes the manufacturing process more complex. However, thanks to the permanent magnets there is no need for rotor flux calculations which makes the mathematical description of the PMSM very simple and also suitable for this type of real-time application on external hardware. The reason for this is that the target hardware, in this case, have much less computing power compared to a desktop computer. Due to the mentioned advantages with the PMSM the focus in this project have been on this motor type to close the loop of model based software development.
Chapter 1. Introduction
1.2 Project Goals
• Deliver a software model of a PMSM suitable for real time execution which also takes magnetic saturation and flux linkage harmonics into consideration.
• Deliver a simple mechanical model for real time execution that simulates external loads and inertia.
• Deliver a fully functional model which is executable on a specific motor drive.
• Deliver a comparison between an inverter using the model and an inverter connected to a real electrical motor.
1.3 System set-up and requirements
The final model should behave as a real electric motor and interface with the already existing control system on a drive. To get a better picture of how the final system should look like at the end, here follows a brief description.
1.3.1 The control system
From a software and control perspective the system of interest can be graphically illus- trated as in Figure 1.1. The controller feeds reference voltages based on the difference be- tween the desired and measured values. The reference signals are then converted from the dq-frame to three phase voltages (urefabs) using the Park transformation and are finally fed into the inverter block. The inverter block consists of all the hardware related to DC/AC inversion.
ref +
Controller dq
Calculations
PMSM
ωm, θm
iabc, uabc
ǫ urefdq
abc
dq abc
urefabc
Inverter
Figure 1.1: Illustration of the drive system controlling an actual motor.
1.4. Thesis outline However, due to the controller being implemented as a vector controller in the dq-frame there was no need to use the final sinusoidal voltages. Instead, the PMSM model could be defined in the dq-frame directly which reduces the number of mathematical operations in the model, i.e no need for Park transformation. In Figure 1.2, the system set-up is shown when replacing the actual motor with the motor model.
ref +
Controller
Calculations
PMSM
ωm, θm
ǫ urefdq
udq, idq
Model
Figure 1.2: Illustration of the drive system controlling a motor model.
1.3.2 Hardware requirements
The target hardware for this application was a STM32F2xx processor running at 200 MHz. However, the loaded software was executed at 4 kHz which means that the time step was limited to 250 ms which is 5 times slower than the recommended time step when performing real-time simulations of electrical systems [1]. Another limitation using this particular CPU was the nonsupport of floating-point data types. This required that the model parameters and variables had to be converted to a fixed-point data type. These issues are explained more in detail in chapter 4.
1.4 Thesis outline
The thesis is separated into six chapters as follows:
• Chapter 1: Presentation of background, purpose and goals.
• Chapter 2: Introducing the theory behind the chosen model approach.
• Chapter 3: Explanation of the finite element analysis done in the project.
• Chapter 4: The Simulink model is explained along with code generation and vali- dation.
Chapter 1. Introduction
• Chapter 5: Presentation of results from the implemented model on the target hard- ware.
• Chapter 6: Discussion and conclusion drawn from the results and suggestions on future work.
Chapter 2 Theory
2.1 Motor equations
The dynamic model of the PMSM used in this project was based on the the well known voltage equations for synchronous machines. The derivation of these equations can be found in many types of publications, e.g., in [2, 3] and the final expressions becomes as follows:
ud= Rsid+ Lddid
dt − ωrLqiq (2.1)
uq = Rsiq+ Lq
diq
dt + ωrLdid+ ωrψP M (2.2) The terminal voltages udq are here expressed by the voltage drop over the stator wind- ings due to the resistance Rs, the voltage drop over the coils which is dependent on the inductance Ldq and the variation in current over time. The last term is a cross-coupling term representing the induced voltage which will appear when the machine is rotating.
In the expression for uq, the fundamental permanent magnet flux ψP M appears due to its condribution to the induced voltage during rotation of the machine. The mentioned set of equations along with a mechanical model could have been sufficient to use in this project which is also shown in [4, 5]. It is also possible to add the saturation effect according to [6]. But considering the application and availability of software, the author chosed an- other method. According to that approach it is convinient to express the terminal voltage in terms of flux-linkage instaed of inductances [7]. This results in the following equations:
ud = Rsid+dψd
dt − ωrψq (2.3)
Chapter 2. Theory
uq = Rsiq+ dψq
dt + ωrψd (2.4)
Using this set of equations together with the relation between flux linkage and current, which is stated in equation 2.5 and 2.6, it is possible to determine the current for the known parameters udq, Ldq, Rsand ψP M.
ψd= ψP M + Ldid (2.5)
ψq = Lqiq (2.6)
When the current and flux linkage are known, the following expression can be used for torque calculation, which also takes the reluctance torque into consideration:
Te = 3p
4 (ψdiq− ψqid) (2.7)
2.2 Complement with magnetic saturation
Why is it necessary to take saturation into acount? The best way to answer that question is to look at the relation between the the total magnetic flux density (B-field) and the ap- plied magnetic field (H-field)(see Figure 2.1). Saturation occurs when the B-field almost stops increasing even with an increased H-field. The magnetic field is dependent on the magnetic properties of the material and will then have an impact on the performance of the machine [7, 8].
0 0.2 0.4 0.6 0.8 1
·104 0
0.5 1 1.5
H [A/m]
B[T]
Figure 2.1: BH-curve for the electrical steel, which the motor in this project was made of [9].
2.3. Adding flux linkage harmonics To consider the magnetic saturation of the motor, equation 2.5 was replaced by FEMM1 generated data for the flux linkage (ψd, ψq) as a function of the current (id, iq).
This means that the current was not calculated in the implemented motor model, instead the current was extracted from a function (see equation 2.8 and 2.9), which was deter- mined using finite element analysis.
id= f (ψd) (2.8)
iq = f (ψq) (2.9)
2.3 Adding flux linkage harmonics
In FEMM is possible to extract the waveform of flux linkage in dq-direction as a function of the rotor position. This can be performed at different current levels but in this project only the no load condition was considered. By performing a spectrum analysis of this function it is possible to extract which orders of harmonics are present in the flux linkage.
By adding this new data to the calculated flux linkage into the model the harmonics have been included in the final model and it can be expressed as follows:
id= f
ψd−
n
X
h=6,12,...
ψd,hcos hpθr 2
(2.10)
iq = f
ψq−
n
X
h=6,12,...
ψq,hsin hpθr 2
(2.11)
As the reader can see, the expressions for the summation of the harmonics (h) are sub- tracted from the fundamental flux linkage. It is also clear that only sixth-order of har- monics are expected which is a result of transforming the three phase flux linkage to dq-quantities.This is shown and further explained in [7].
2.4 Mechanical model
The mechanical part of the motor was modelled using a common and simple model which is shown in equation 2.12. By solving this equation for ωm for a certain inertia J, load torque Tl, friction and produced torque Te it is possible to determine the corresponding speed.
1Finite Element Method Magnetics is an open source simulation software developed by Dr. David Meeker :http://www.femm.info/wiki/HomePage
Chapter 2. Theory
Te− Tl = Jdωm
dt + f riction (2.12)
The friction was modelled using equation 2.13 and the reason was that it would then be possible to have a constant friction together with a speed dependent friction. It is then possible to tune the friction to make it more realistic and it could also be a complement to the constant load torque acting as a load.
f riction = ωmb + sign(ωm)c (2.13) Note that the notation for mechanical rotor speed is ωm and electrical rotor speed is ωr
where the electrical rotor speed is used in the already stated motor equations.
ωr = p
2ωm (2.14)
Chapter 3
Finite Element Method Analysis
The finite element analysis (FEA) in his project was carried out using a software called Finite Element Method Magnetics (FEMM). This is an open source software devel- oped by Dr. David Meeker and is suitable for low frequency 2D FEA [10]. In this chapter all steps from geometry design to data extraction are explained.
3.1 The geometry
To run simulations and extract usable data from FEMM, the first thing to do was to make a geometry of the target object. In this project the target was based on drawings of a 10kW PMSM that was available for load tests and back-emf measurements. The specific geometrical numbers will not be stated in this report but in Table 3.1 some the used motor parameters are presented and in Figure 3.1 the final FEMM model is shown.
Number of stator slots (Qs) 18
Number of poles (p) 8
Rs 17.9 mOhm
Ld 0.38 mH
Lq 0.79 mH
Table 3.1: Motor parameters from data sheets and drawings
As shown in Figure 3.1, the stator winding is divided into two layers meaning that the coils are pitched. Each phase consists of six coils in total where three are connected in series (see Figure 3.2). In FEMM this was modelled by divided each phase into two circuits with three coils in series. The current running through the circuit was then defined as half the phase current.
Chapter 3. Finite Element Method Analysis
N S
Figure 3.1: The geometry of the PMSM in FEMM with winding arrangement where brown, black and gray corresponds to phase A,B and C.
A
B C
Figure 3.2: Illustration of the Y connected coils
3.2 Extrapolation of BH-curve
The magnetic properties of the steel laminations have a big impact of when investigating magnetic saturation in the machine. In most cases, there are only measurements up to a level of 1.8 T in flux density but the machine can be operating at higher levels. The BH- curve was also extrapolated a bit to give FEMM a few more data points at higher field intensities. This was done using Law of Approach to Saturation Extrapolation (LAS) which according to [11] is a sufficient method. The idea of LAS is to solve a system of equations (see equations 3.1) for known B, H and µ0 and solve for the saturation
3.3. Validation magnetization (Ms) and the curve fitting coefficient b. By substituting the calculated Ms
and b into 3.1 provides an analytical expression for the saturation of the B-field.
B− µ0
H + Ms− MsHb2
dB
dH = µ0
1 + 2MsHh3
(3.1)
In Table 3.2, the resulting Ms and b are shown together with the last data points for the magnetic filed used in equation 3.1. The extrapolation was done up to 2.2 T which was considered high enough for further simulations. In Figure 3.3, the extrapolated part is plotted together with the original part of the BH-curve.
Blast 1.8 [T]
Hlast 9390 [A/m]
µ0 4π· 10−7 [H/m]
Ms 1.4319· 106 [A/m]
b 5.4481· 106 [A2/m2]
Table 3.2: Blast and Hlast is the last data points from data sheet which resulted in the calculated parameters Ms and b
0 0.5 1 1.5
·105 0
0.5 1 1.5 2
H [A/m]
B[T]
Extrapolated Original
Figure 3.3: Extrapolated BH-curve for the magnetic steel in the machine up to a magnetic field of 2.0 T
3.3 Validation
To validate the FEMM model, the theory of electromotive force (EMF) was used. Due to the constantly active flux from the permanent magnets, it is possible to measure the back- EMF by rotating the machine and measure the terminal voltage. In the same way, it is
Chapter 3. Finite Element Method Analysis
possible to calculate the resulting flux for each phase for different rotor angles in FEMM.
This calculation resulted in the function Ψ(θm). To transform this function to voltage, Faraday’s law of induction was used which says that:
= dΨ
dt (3.2)
In this case, the flux was expressed as a function of the rotor angle which made it conve- nient to express the time derivative as following.
= dΨ dθm · dθm
dt (3.3)
Hence, the back-EMF of the FEMM model was calculated by determining the rotor angle derivative of the flux and multiply it with the desired speed. The measured back-EMF from the real machine was taken by connecting an oscilloscope to the terminal connectors and then spin the rotor which resulted in an induced voltage on the terminal side. The comparison between FEMM calculated and measured back-EMF is shown in Figure 3.4 and corresponds to a speed of 2000 RPM.
0 20 40 60 80
−200 0 200
Mechanical degrees θm[Deg]
Line-to-linevoltageUAB[V] FEMM Measured
Figure 3.4: Comparison between FEMM calculated and measured back-EMF at 2000 rpm
3.4 Data extraction
As described in chapter 2 the idea with the chosen model approach was to also take magnetic saturation into consideration. In other words, it was necessary to investigate how the flux linkage was related to the applied current. In FEMM it was possible to extract data of the flux linkage as function of current (ψ(i)) by run simulations and vary id and iq respectively. The current was varied between 500 A and -500 A and the results of these simulations are shown in Figure 3.5a and 3.5b.
3.4. Data extraction
−600 −400 −200 0 200 400 600 0
0.1 0.2
id[A]
ψd[Vs]
(a)
−600 −400 −200 0 200 400 600
−0.2
−0.1 0 0.1 0.2
iq[A]
ψq[Vs]
(b)
Figure 3.5: The extracted flux linkage where, in Figure 3.5a, −500 ≤ Id ≤ 500, Iq = 0 and, in Figure 3.5b, −500 ≤ Iq ≤ 500, Id= 0.
As described in chapter 2, the magnitude of the flux linkage harmonics was also considered. By extracting the no-load flux linkage in d-,q-direction for 0 ≤ θm ≤ 2π it was possible to see the behavior of ψd and ψq for different rotor angles. The results of these simulations are shown in Figure 3.6.
0 0.5 1 1.5
0.1686 0.1688 0.169 0.1692 0.1694 0.1696 0.1698
θr[rad]
ψd[Vs]
(a)
0 0.5 1 1.5
0 0.2 0.4 0.6 0.8 1 1.2
·10−3
θr[rad]
ψq[Vs]
(b)
Figure 3.6: No-load flux linkage in d and q direction for 0 ≤ θm ≤ 2π
The extracted data shown in Figure 3.6 was then analysed using Fourier transfor- mation. This was carried out using existing functions in Mathlab and the results are shown in Figure 3.7. As described in chapter 2, only the 6th order of harmonics was expected.
To verify this, a fundamental frequency of 50 hz was chosen and as the reader can see, the peaks appear at frequencies which are 6, 12 and 18 times the fundamental frequency.
Chapter 3. Finite Element Method Analysis
200 400 600 800 1,000 0
2 4 6 ·10−4
Frequency [Hz]
Amplitude[Vs]
(a)
200 400 600 800 1,000 0
1 2 3 4
·10−4
Frequency [Hz]
Amplitude[Vs]
(b) Figure 3.7: Frequency spectrum of no load flux linkage
Chapter 4 Simulink
To simulate the behavior of the presented mathematical motor model and also to convert it to C-code, a Mathworks software called Simulink was used. This software makes it possible to solve the equations using block-programming where the blocks are predefined mathematical operations, e.g., integrals, multiplication, addition, etc. This soft- ware makes it also possible to convert the model into code which then can be used as an executable function one external hardware. This loop from model to code are commonly called model based software development and this chapter consist of the steps which were taken in this project.
4.1 The continuous model
The first step was to create a reference model to ensure that the simulations converged to an expected value. This model was based on the assumption that the simulations were im- plemented in continuous time. Here, both the linear and the saturated flux linkage model was implemented and it was possible to switch between them for comparison. As the reader can see in Figure 4.1, 4.2 and 4.3, the equations from chapter 2 was simply de- fined in Simulink and reference simulations was done with provided motor parameters (see Table 3.1). In Simulink, it is also possible to explicit define which solver to use. In this stage that option was set to ”automatic” which allows Simulink to choose the most suitable solver for the particular problem. By calculating steady-state values, i.e., the dif- ferential terms in the equations set to 0, and compare them to the steady-state values from the simulations, it was found that the model performed as expected.
Chapter 4. Simulink
Figure 4.1: Current (id) calculations
Figure 4.2: Calculations of the produced torque (Te).
4.1. The continuous model
Figure 4.3: Continuous mechanical load model for rotor speed and position calculations based on produced torque, load torque, inertia and friction. The switch makes it possible to bypass the load model and set a constant speed.
To include the saturation effect in the model, the FEM-generated data of the flux linkage was stored in look-up tables and was also inverted from i(ψ(I)) to I(ψ). However, it will later be discovered that this inversion of the flux-linkage function had a big impact on the resulting currents and torque. During the first simulations with magnetic saturation the saturation curve was simply inverted by change place of Y- and X-axis (see Figure 4.4). This resulted in an even higher produced torque compared to just using the provided motor parameters, i.e., the linear model. To investigate this even further the current was limited to a maximum value which shown that the flux linkage could reach high levels where the corresponding current then would be even greater compared to the linear model.
Figure 4.4: The simple inversion methods of the saturation curve
Chapter 4. Simulink
By introducing another invertion method found in [12] this problem was solved (see Figure 4.5). In practice there is no bigger difference between the two inversion methods but due to the algebraic constraint in the second method, Simulink is forced to take small steps when extracting the data from the look-up table. This resulted in a distinct limitation of the flux linkage and also the produced torque which can be seen in Figure 4.6. In the same module of the model, the flux linkage harmonics were also introduced and added to the calculated flux linkage based on the rotor position (see Figure 4.5). The resulting rippled current along with the corresponding smooth current can be seen in Figure 4.7.
When analysing the resulting torque in Figure 4.8, the harmonic components in iqare also clearly appearing in the torque. This is the result of the product of iqand ψdin the torque equation which shows that current in the q-direction is highly affecting the produced torque.
Figure 4.5: The advanced inversion methods of the saturation curve
4.1. The continuous model
0 0.5 1 1.5 2 2.5 3 0
100 200 300
Time [s]
Torque[Nm]
Saturated Unsaturated
(a)
0 0.5 1 1.5 2 2.5 3 0
5· 10−2 0.1 0.15 0.2
Time [s]
ψ[Vs]
ψq,Saturated
ψd,Saturated
ψd
ψq
(b)
Figure 4.6: Locked rotor simulation of torque (a) and flux linkage (b) with and without saturation at ˆU = 5V and θu = 135◦
2.8 2.85 2.9 2.95 3
−12
−11
−10
−9
Time [s]
Id[A]
(a)
2.8 2.85 2.9 2.95 3 7
7.2 7.4
Time [s]
Iq[A]
(b)
Figure 4.7: Simulation of Id(a) and Iq(b) with ripple due to flux linkage harmonics
Chapter 4. Simulink
2.8 2.85 2.9 2.95 3
7.3 7.4 7.5 7.6 7.7
Time [s]
Torque[Nm]
Figure 4.8: Simulation of torque with ripple due to flux linkage harmonics
4.2 Discretization methods
The final software model will be run as a function which will be called every 250 µs, i.e., 4 kHz. Hence, the input and output signals will be updated with a fixed time step, the definition of time in the model needed to be changed from continuous to discrete. In Simulink, there are three different discretization methods to choose from: forward Euler, backward Euler and trapezoidal discretization. The simplest of these methods is the for- ward Euler method which is an explicit method capable of generating accurate solutions using fairly low computing power. However, the most crucial drawback with this method is the fact that it is highly dependent on the step length. Solving differential equation using this method requires a very small time step to provide a numerical stable solution [13].
This method was implemented and tested and it was found that for even very short time steps (down to 25µs), it was impossible to perform stable simulations of the model. To overcome this problem the method was changed to the backward Euler method. This im- plicit method can generate numerical stable solutions for larger time steps compared to the forward Euler method. But every improvement comes with a price which in this case affected the needed computing power and the signal routing through the model. Using this method, the integration block does no longer brake algebraic loops which then had to be considered. To solve this, delays were added into the model to brake the algebraic loops.
The model was then simulated and a stable solution was found at time steps around 50 µs which was still 5 times shorter than the required 250 µs. When integrated in the software platform, this would require that the model must be executed at least 5 times for every 250 µs cycle. The other problem using this method is the more advanced implementation when solving implicit equations (backward Euler) compared to solving explicit equations
4.3. Fixed point conversion (forward Euler) [13, 14]. However, to conclude if the model would be too heavy for the target system the actual code needed to be generated and tested on the target hardware.
4.3 Fixed point conversion
To generate C-code from a Simulink model for real-time execution, all signals in the model needed to be converted from the default data type, double, to a fixed-point data type. What this means is that variables and constants should be represented by a signed or unsigned integer with a number of bits where a scaling factor determines how much one bit represents. This does also mean that if using an insufficient scaling factor on a variable, the integer can overflow or generate resolution problems resulting in unexpected behaviors in the model. The fixed point conversion can be done using a Simulink tool called Fixed-Point Tool where the system, based on simulations, proposes a scale factor for a given bit length. However, in the software platform, a scaling for certain units was already defined and to simplify the interfacing between the motor model and the control system, the existing scaling for 16 bit integers was used. Through simulation, the conver- sion to fixed-point data types was analyzed to ensure that the variables did not overflow or had lack of resolution leading to big rounding errors.
4.4 Code generation and validation
The C-code generation was carried out using the Simulink add-on package Embedded Coder and by following the manual for code generation [15]. Based on hardware tar- get specific parameters, Embedded Coder could generate code suitable for the hardware used in this project with input and output data interfaces that matched the target software platform. The code validation was done in two steps. The first step was to compile and execute the code in Simulink using a procedure called software in the loop (SIL) [16].
In practice, the code is generated as a Simulink block (see Figure 4.9) which makes it possible to run the code and compare it to its original. In Figure 4.10 one result of the comparison between the model and SIL block can be seen. As expected, the error is zero which reveals that the generated code has not been optimized in a way that affects the result. The second step in the process of validation was to execute the generated code on the target hardware where the produced torque was compared with FEMM simulations and also measurements of a real machine. These results are presented in chapter 5.
Chapter 4. Simulink
Figure 4.9: A SIL block consisting of generated C-code
4.4. Code generation and validation
Figure 4.10: Comparison of iq between the discrete motor model and the generated C code
Chapter 4. Simulink
Chapter 5 Results
The implementation of the generated code resulted in a real-time executable model on the target hardware interfaced with the on board control system. The performance analysis was made by comparing FEMM simulations, the software model, Table 5.1, and the back- to-back mounted PMSM, Table 5.2, in terms of produced torque. The comparison was carried out by first taking measurements from the actual PMSM where the currents id
and iq was extracted for a number of fixed values of the produced torque. The same currents were then used as input in FEMM simulations where the corresponding torque was calculated using Maxwell stress tensor in the middle of the air gap for different rotor positions. In Table 5.1, the reader can see both the measured value and the mean value for different rotor positions of the simulated torque extracted from FEMM and also the difference in percent.
id[A] iq[A] Tmeasured[Nm] TF EM [Nm] Diff [%]
-44 114 100 127.0 ∼27
-99 205 200 238.7 ∼19
-160 280 300 332.9 ∼11
Table 5.1: Comparison between measured and FEM simulated torque
To investigate the performance of the implemented software model, the same cur- rents was set as reference currents in the control system on the inverter. In Table 5.2 is the FEMM simulated torque compared to the reported torque from the real time software model.
Chapter 5. Results
id[A] iq[A] TF EM [Nm] Treal−timemodel[Nm] Diff [%]
-44 114 127.0 126.0 <1
-99 205 238.7 225 ∼6
-160 280 332.9 302.6 ∼10
Table 5.2: Comparison between measured and FEM simulated torque
As a complement to the earlier mentioned torque comparison, some other perfor- mance tests were done by applying load steps and speed steps in the controller. Note that these tests were performed to see that the software model was interfacing the control sys- tem in a sufficient way rather than a performance test of the actual controller. Hence, the controller was not optimized for this typical motor, just roughly tuned so it could reach a reference speed. The first test was done by applying a load step to the load torque variable in the software model. The reference speed was set to 2000 rpm at no load and a load step of 150 Nm was then performed in the mechanical part of the model (see Figure 5.1).
0 0.2 0.4 0.6 0.8 1 0
500 1,000 1,500 2,000
Time [s]
Amplitude[rpm] Actual speed Reference Speed
Load torque Actual torque
(a)
0 0.2 0.4 0.6 0.8 1 0
1,000 2,000
Time [s]
Amplitude[rpm]
Actual speed Reference Speed
Load torque Actual total current
(b)
Figure 5.1: Results of load step test in speed (a) and current (b). Note that the scale on y-axis is in rpm and the current in (b) an unscaled 16 bit signed integer
Chapter 6
Discussion and Conclusions
6.1 Discussion
When analyzing the modelling approach and the results, the outcome can be seen from two perspectives. From an application perspective where the goal was to develop a real time software model which should behave as a PMSM and generate necessary feedback to the control system, the results are good enough. However, from an electrical machine analysis perspective, the results may leave a little bit more to wish for. To explain what this mean, it is convenient to divide the project into two parts, FEMM and Simulink.
6.1.1 FEMM
When comparing the FEMM model with the real PMSM in terms of produced torque, the FEMM model does produce more torque for the same amount of current as the actual machine. The answer to that lay probably in the geometry of the FEMM model. The presence of air in the machine have a big impact on the produced torque and there are a lot of production tolerances which could make a differences in the result. This was considered by adding 0.3 mm of air below the permanent magnets which improved the result. An actual machine has also losses due to the stacking and welding of the thin steel laminates which are not considered in the FEA due to the magnetic solutions only were defined in two dimensions. The last factor that probably also had a big impact was that the cross saturation was not considered when extracting saturation data from the model. This procedure was considered as too time consuming and was then neglected. But according to specialists at Inmotion, the cross saturation has a quite big impact on this particular motor type.
Chapter 6. Discussion and Conclusions
6.1.2 Simulink
A big part of the final implemented software model was based on data extracted from FEMM. Hence, the real time model was expected to be very close to the FEMM model in terms of produced torque. According to the results in table 5.2 this was also the case even if the error is growing for higher levels of current. However, the error stays more or less below 10 % which is within the expected range. But it would of course be interesting to have more data points to investigate the differences for even higher levels of current.
A guess is the error will grow even more and the reason is probably how the FEMM data was integrated in the simulink model. As explained in chapter 4 the process of inverting the function ψ(i) into i(ψ) had a big impact on the result. Due to usage of unsupported blocks it was not possible to use the more advanced inversion method from the continuous model when converting the model to a discretized, fixed-point model.
6.2 Conclusions
In this thesis the modelling process to generate a real-time software model of an actual PMSM using a modelling approach including FEA and Simulink simulations has been presented. The code generated real-time model has been sufficiently executed on the tar- get hardware interfacing the on board control system. The model was also taking magnetic saturation into consideration but flux linkage harmonics had to be neglected due to sta- bility issues in the final discrete model. A comparison of the performance of the software model and the real PMSM have been done in terms of produced torque. This comparison resulted in that the modelling process need refinements to come closer to reality even if the stated results were considered as good enough for this application.
6.3 Suggestions on future work
6.3.1 Improvement of the FEMM-model
To improve the FEA of the PMSM, a 3D model of the motor could be implemented. This would also take end effects and a third magnetic field direction into account. This model could then be used to generate more accurate inductance or flux linkage tables to improve torque estimations in the control system.
6.3.2 Develop the mechanical model
The mechanical model implemented in this project have some degrees of variability where both inertia and load torque can be varied. It would be useful for users to be able to insert data from the actual vehicle environment to make performance tests of the tuned control
6.3. Suggestions on future work system. One example of this could be the possibility to insert driving cycles of a forklift where also the carried load is variable.
6.3.3 Modelling of other motor types (IM)
Due to the fact that Inmotions products also are used for controlling induction machines, it would also be very useful to have a software model for this type of electrical machine as well. It can be noted that this was, from the beginning, one part of this project which resulted in one FEMM model and a Simulink model. Due to some major problems regard- ing stability in the discrete Simulink model and also verification of the FEMM model, this part of the project was closed and left as a future project.
6.3.4 Study the fixed-point conversion work flow
As mentioned in the section describing the fixed-point conversion of the Simulink model, there was possible to use Simulink’s own tool for this conversion of data types. However, the authors impression, after using this tool as one part of the model based software de- velopment, was that the tool is closely related to a black box. All signals are automatically scaled based the defined size or fraction length of the signal to fit a defined range. This re- sults in many unique scale factors which can be a problem if the signal should be used or displayed in other parts of the system. This is why this tool must be further investigated if it is going to be used professionally together with the model based software development approach.
Chapter 6. Discussion and Conclusions
Appendix A. Simulink models
Appendix A
Simulink models
A.0.1 Continuous Model
32
Appendix A. Simulink models
Appendix A. Simulink models
A.0.2 Discrete Model
34
Appendix B
Generated C code
B.0.1 Main file
/*
* Academic License - for use in teaching, academic research, and meeting
* course requirements at degree granting institutions only. Not for
* government, commercial, or other organizational use.
*
* File: PM_motor0.c
*
* Code generated for Simulink model ’PM_motor0’.
*
* Model version : 1.42
* Simulink Coder version : 8.9 (R2015b) 13-Aug-2015
* C/C++ source code generated on : Wed Nov 16 13:31:06 2016
*
* Target selection: ert.tlc
* Embedded hardware selection: STMicroelectronics->ST10/Super10
* Code generation objectives:
* 1. Execution efficiency
* 2. RAM efficiency
* Validation result: Not run
*/
#include "PM_motor0.h"
#ifndef UCHAR_MAX
#include <limits.h>
#endif
#if ( UCHAR_MAX != (0xFFU) ) || ( SCHAR_MAX != (0x7F) )
Appendix B. Generated C code
#error Code was generated for compiler with different sized uchar/
char. \
Consider adjusting Test hardware word size settings on the \ Hardware Implementation pane to match your compiler word sizes as
\
defined in limits.h of the compiler. Alternatively, you can \ select the Test hardware is the same as production hardware option
and \
select the Enable portable word sizes option on the Code Generation > \
Verification pane for ERT based targets, which will disable the \ preprocessor word size checks.
#endif
#if ( USHRT_MAX != (0xFFFFU) ) || ( SHRT_MAX != (0x7FFF) )
#error Code was generated for compiler with different sized ushort /short. \
Consider adjusting Test hardware word size settings on the \ Hardware Implementation pane to match your compiler word sizes as
\
defined in limits.h of the compiler. Alternatively, you can \ select the Test hardware is the same as production hardware option
and \
select the Enable portable word sizes option on the Code Generation > \
Verification pane for ERT based targets, which will disable the \ preprocessor word size checks.
#endif
#if ( UINT_MAX != (0xFFFFU) ) || ( INT_MAX != (0x7FFF) )
#error Code was generated for compiler with different sized uint/
int. \
Consider adjusting Test hardware word size settings on the \ Hardware Implementation pane to match your compiler word sizes as
\
defined in limits.h of the compiler. Alternatively, you can \ select the Test hardware is the same as production hardware option
and \
select the Enable portable word sizes option on the Code Generation > \
Verification pane for ERT based targets, which will disable the \ preprocessor word size checks.
#endif
#if ( ULONG_MAX != (0xFFFFFFFFUL) ) || ( LONG_MAX != (0x7FFFFFFFL) )
#error Code was generated for compiler with different sized ulong/
long. \
Appendix B. Generated C code
Consider adjusting Test hardware word size settings on the \ Hardware Implementation pane to match your compiler word sizes as
\
defined in limits.h of the compiler. Alternatively, you can \ select the Test hardware is the same as production hardware option
and \
select the Enable portable word sizes option on the Code Generation > \
Verification pane for ERT based targets, which will disable the \ preprocessor word size checks.
#endif
/* Block signals and states (auto storage) */
DW rtDW;
/* External inputs (root inport signals with auto storage) */
ExtU rtU;
/* External outputs (root outports fed by signals with auto storage) */
ExtY rtY;
/* Real-time model */
RT_MODEL rtM_;
RT_MODEL *const rtM = &rtM_;
extern int16_T look1_is16lu32n16ts16D_gNlSsjED(int16_T u0, const int16_T bp0[],
const int16_T table[], uint32_T maxIndex);
extern void mul_wide_s32(int32_T in0, int32_T in1, uint32_T * ptrOutBitsHi,
uint32_T *ptrOutBitsLo);
extern int32_T mul_s32_s32_s32_sr15_sat(int32_T a, int32_T b);
extern int32_T mul_s32_s32_s32_sr46(int32_T a, int32_T b);
extern int32_T mul_s32_s32_s32_sr30_sat(int32_T a, int32_T b);
extern int32_T div_nzp_repeat_s32_sat_floor(int32_T numerator, int32_T
denominator, uint16_T nRepeatSub);
extern uint32_T div_nzp_repeat_u32_sat_ceiling(uint32_T numerator, uint32_T
denominator, uint16_T nRepeatSub);
extern uint32_T div_nzp_repeat_u32_sat(uint32_T numerator, uint32_T denominator,
uint16_T nRepeatSub);
int16_T look1_is16lu32n16ts16D_gNlSsjED(int16_T u0, const int16_T bp0[], const
int16_T table[], uint32_T maxIndex) {
uint32_T frac;
Appendix B. Generated C code
uint32_T iRght;
uint32_T iLeft;
/* Lookup 1-D
Canonical function name: look1_is16lu32n16ts16Ds32_binlcs Search method: ’binary’
Use previous index: ’off’
Interpolation method: ’Linear’
Extrapolation method: ’Clip’
Use last breakpoint for index at or above upper limit: ’off’
Remove protection against out-of-range input in generated code: ’off’
Rounding mode: ’simplest’
*/
/* Prelookup - Index and Fraction Index Search method: ’binary’
Extrapolation method: ’Clip’
Use previous index: ’off’
Use last breakpoint for index at or above upper limit: ’off’
Remove protection against out-of-range input in generated code: ’off’
Rounding mode: ’simplest’
*/
if (u0 <= bp0[0UL]) { iLeft = 0UL;
frac = 0UL;
} else if (u0 < bp0[maxIndex]) { /* Binary Search */
frac = maxIndex >> 1UL;
iLeft = 0UL;
iRght = maxIndex;
while (iRght - iLeft > 1UL) { if (u0 < bp0[frac]) {
iRght = frac;
} else {
iLeft = frac;
}
frac = (iRght + iLeft) >> 1UL;
}
frac = ((uint32_T)((uint16_T)u0 - bp0[iLeft]) << 16) / ((
uint16_T)bp0[iLeft + 1UL] - bp0[iLeft]);
} else {
iLeft = maxIndex - 1UL;
frac = 65536UL;
}
Appendix B. Generated C code
/* Interpolation 1-D
Interpolation method: ’Linear’
Use last breakpoint for index at or above upper limit: ’off’
Rounding mode: ’simplest’
Overflow mode: ’wrapping’
*/
return (int16_T)(((int32_T)table[iLeft + 1UL] - table[iLeft]) * (int32_T)frac >>
16) + table[iLeft];
}
void mul_wide_s32(int32_T in0, int32_T in1, uint32_T *ptrOutBitsHi , uint32_T
*ptrOutBitsLo) {
uint32_T absIn0;
uint32_T absIn1;
uint32_T in0Lo;
uint32_T in0Hi;
uint32_T in1Hi;
uint32_T productHiLo;
uint32_T productLoHi;
absIn0 = (uint32_T)(in0 < 0L ? -in0 : in0);
absIn1 = (uint32_T)(in1 < 0L ? -in1 : in1);
in0Hi = absIn0 >> 16UL;
in0Lo = absIn0 & 65535UL;
in1Hi = absIn1 >> 16UL;
absIn0 = absIn1 & 65535UL;
productHiLo = in0Hi * absIn0;
productLoHi = in0Lo * in1Hi;
absIn0 *= in0Lo;
absIn1 = 0UL;
in0Lo = (productLoHi << 16UL) + absIn0;
if (in0Lo < absIn0) { absIn1 = 1UL;
}
absIn0 = in0Lo;
in0Lo += productHiLo << 16UL;
if (in0Lo < absIn0) { absIn1++;
}
absIn0 = (((productLoHi >> 16UL) + (productHiLo >> 16UL)) + in0Hi * in1Hi) +
absIn1;
Appendix B. Generated C code
if (!((in0 == 0L) || ((in1 == 0L) || ((in0 > 0L) == (in1 > 0L))) )) {
absIn0 = ˜absIn0;
in0Lo = ˜in0Lo;
in0Lo++;
if (in0Lo == 0UL) { absIn0++;
} }
*ptrOutBitsHi = absIn0;
*ptrOutBitsLo = in0Lo;
}
int32_T mul_s32_s32_s32_sr15_sat(int32_T a, int32_T b) {
int32_T result;
uint32_T u32_chi;
uint32_T u32_clo;
mul_wide_s32(a, b, &u32_chi, &u32_clo);
u32_clo = u32_chi << 17UL | u32_clo >> 15UL;
u32_chi = (uint32_T)((int32_T)u32_chi >> 15L);
if (((int32_T)u32_chi > 0L) || ((u32_chi == 0UL) && (u32_clo >=
2147483648UL))) {
result = MAX_int32_T;
} else if (((int32_T)u32_chi < -1L) || (((int32_T)u32_chi == -1L ) && (u32_clo <
2147483648UL))) { result = MIN_int32_T;
} else {
result = (int32_T)u32_clo;
}
return result;
}
int32_T mul_s32_s32_s32_sr46(int32_T a, int32_T b) {
uint32_T u32_chi;
uint32_T u32_clo;
mul_wide_s32(a, b, &u32_chi, &u32_clo);
return (int32_T)u32_chi >> 14L;
}
int32_T mul_s32_s32_s32_sr30_sat(int32_T a, int32_T b) {
int32_T result;
Appendix B. Generated C code
uint32_T u32_chi;
uint32_T u32_clo;
mul_wide_s32(a, b, &u32_chi, &u32_clo);
u32_clo = u32_chi << 2UL | u32_clo >> 30UL;
u32_chi = (uint32_T)((int32_T)u32_chi >> 30L);
if (((int32_T)u32_chi > 0L) || ((u32_chi == 0UL) && (u32_clo >=
2147483648UL))) {
result = MAX_int32_T;
} else if (((int32_T)u32_chi < -1L) || (((int32_T)u32_chi == -1L ) && (u32_clo <
2147483648UL))) { result = MIN_int32_T;
} else {
result = (int32_T)u32_clo;
}
return result;
}
int32_T div_nzp_repeat_s32_sat_floor(int32_T numerator, int32_T denominator,
uint16_T nRepeatSub) {
int32_T quotient;
uint32_T absNumerator;
uint32_T absDenominator;
absNumerator = (uint32_T)(numerator >= 0L ? numerator : - numerator);
absDenominator = (uint32_T)(denominator >= 0L ? denominator : - denominator);
if ((numerator < 0L) != (denominator < 0L)) {
absNumerator = div_nzp_repeat_u32_sat_ceiling(absNumerator, absDenominator,
nRepeatSub);
quotient = absNumerator <= 2147483647UL ? -(int32_T) absNumerator :
MIN_int32_T;
} else {
absNumerator = div_nzp_repeat_u32_sat(absNumerator, absDenominator,
nRepeatSub);
quotient = absNumerator <= 2147483647UL ? (int32_T) absNumerator :
MAX_int32_T;
}
return quotient;
Appendix B. Generated C code
}
uint32_T div_nzp_repeat_u32_sat_ceiling(uint32_T numerator, uint32_T denominator,
uint16_T nRepeatSub) {
uint32_T quotient;
uint16_T iRepeatSub;
boolean_T numeratorExtraBit;
boolean_T overflow;
overflow = false;
quotient = numerator / denominator;
numerator %= denominator;
for (iRepeatSub = 0U; iRepeatSub < nRepeatSub; iRepeatSub++) { numeratorExtraBit = (numerator >= 2147483648UL);
numerator <<= 1UL;
overflow = (overflow || (quotient >= 2147483648UL));
quotient <<= 1UL;
if (numeratorExtraBit || (numerator >= denominator)) { quotient++;
numerator -= denominator;
} }
if (overflow || (quotient >= MAX_uint32_T)) { quotient = MAX_uint32_T;
} else {
if (numerator > 0UL) { quotient++;
} }
return quotient;
}
uint32_T div_nzp_repeat_u32_sat(uint32_T numerator, uint32_T denominator,
uint16_T nRepeatSub) {
uint32_T quotient;
uint16_T iRepeatSub;
boolean_T numeratorExtraBit;
boolean_T overflow;
overflow = false;
quotient = numerator / denominator;
numerator %= denominator;
for (iRepeatSub = 0U; iRepeatSub < nRepeatSub; iRepeatSub++) { numeratorExtraBit = (numerator >= 2147483648UL);
