IPM Motor Simulation and FOC Software

[Pete9008]

With the number of problems people seem to be having with the field weakening in the FOC software I'm getting a bit uncomfortable spending a lot of time on a Prius Gen3 inverter system install and am wondering whether to move to the Outlander inverter (I'm using the Outlander rear motor). I'd prefer not to though as the OpenInverter/Gen3 Inverter combination would be a very nice solution. To try and figure out which way to go I've decided to try to put together a crude motor model and simulation program to try and better understand what is going on.

So far I have a basic motor model running in a Qt program that graphs the results. The model just uses the Park/Clarke transforms to go from the inverter PWM voltages to Vq and Vd, uses these along with a calculated BEMF and Ld and Lq to calculate motor currents, uses the currents to calculate torque, the torque to calculate acceleration, the acceleration to calculate the change in rotor position and then the inverse Park/Clarke to get back to 3phase motor currents. No idea whether this is a valid model (and anything more sophisticated involving magnetic field modelling is well beyond me) but so far it seems to generate vaguely sensible results.

The next step is to add the OpenInverter FOC software and see if the model continues to give sensible results and whether the field weakening issues can be reproduced.

So the question - what's the bare minimum subset of the OpenInverter FOC software that I need to integrate into the simulation? It needs to be enough to be realistic but not so much that I spend forever writing test stubs for code I don't really understand.

EDIT (Jan 23)
This thread has become quite a lengthy tome. If you are interested in reading through the full process of my getting to grips with how motors work and lots of plots please feel free to read on. However if you are more interested in a high level understanding and using the simulator without having to read through the whole thread here are links to the posts that are hopefully of the most interest

Wiki link - https://openinverter.org/wiki/IPM_Motor_Simulator

IPM Equations - viewtopic.php?p=44359#p44359
How OpenInverter does Space Vector Modulation - viewtopic.php?p=44564#p44564
MTPA curves - viewtopic.php?p=45061#p45061
Explanation of SyncAdv change - viewtopic.php?p=45088#p45088
Explanation of PM Motors - viewtopic.php?p=45449#p45449
Explanation of control loops/tuning - viewtopic.php?p=45510#p45510
Start of discussion about motor inductance - viewtopic.php?p=46915#p46915
Motor equations voltage terms - https://openinverter.org/forum/vie ... 2#p47212
Start of results using FULL motor model - viewtopic.php?p=47226#p47226 (results prior to this point can be misleading!)
Current circles and voltage ellipses - viewtopic.php?p=47287#p47287
First coasting test for critical current - viewtopic.php?p=47564#p47564
Power/Torque plots - viewtopic.php?p=47770#p47770
Comparison of simulator vs real world - viewtopic.php?p=48144#p48144
Johannes says 'topic solved' - viewtopic.php?p=48417#p48417
Calculating flux linkage from plots - viewtopic.php?p=48609#p48609
Over modulation - viewtopic.php?p=48638#p48638
Thoughts for future development - viewtopic.php?p=48756#p48756
Johannes video on how to use simulator - viewtopic.php?p=48960#p48960
Vq limiting (not resolved) - viewtopic.php?p=49172#p49172
How to install - viewtopic.php?p=49176#p49176 viewtopic.php?p=50233#p50233 viewtopic.php?p=50232#p50232
Full install guide - viewtopic.php?p=50884#p50884
Operating Instructions - viewtopic.php?p=50885#p50885

[johu]

Great effort! I think you just need the code in PwmGeneration::Run() and all the functions it calls. And probably PwmGeneration::SetTorquePercent() as that does the MTPA calculation. Is the model viewable anywhere yet?

[Pete9008]

Thanks, that's great. Is it just the PwmGeneration class that's needed, does it also include the PID loops for Iq, Id and FW?

I'll have a go at adding it this afternoon.

Once I have a little more confidence in the code it I'll put it up. Any suggestion on how/where (I don't get on well with GIT and don't have an account!).

[Pete9008]

Also which version of the code is it best to start with?

[Bigpie]

Looking forward to taking a look at this simulation. Github would be great, but failing that Google Drive or something?

I'm currently trying to get my head around MPTA and measuring motor electrical parameters to replace the hardcoded values to see if torque improves appreciably before attempting to make them editable from the web interface.

viewtopic.php?t=2598

[RetroZero]

Thanks for having a stab at this problem. Way over my pay grade

[johu]

Also which version of the code is it best to start with?

Most recent master branch
https://github.com/jsphuebner/stm32-sin ... on-foc.cpp

You will see dependencies to other modules there, most importantly foc.cpp and picontroller.cpp which reside in libopeninv

[Pete9008]

Thanks, not got as far as I'd hoped this afternoon and haven't managed to look at integrating the OpenInverter code yet. Fixed a couple of bugs in the model but mainly been fighting with Qt graphs; they were taking over a minute to render a seconds worth of simulation! Eventually found the problem was that it was re-rendering the graph as each point was added - changed it to only render complete blocks and it's fine now (renders 10sec worth of data in a couple of seconds - may become a problem again as more lines are added though).

This is where it is currently at:

The control windows allows the model parameters to be set and control running the simulation. It currently pops up two windows, one for motor current and one for motor speed. This simulation is just starting with a motor spinning at 3600rpm feeding into a 1R load and watching how it spins down - it looks fairly sensible. Changing the parameters and load resistance change things in sensible ways which is promising.

Currently it doesn't do anything about core saturation (the simulation gives 18000A into a short circuit!), the drag parameter doesn't do anything and I've just realised that I've forgotten to take account of the number of motor poles.

Also bear in mind that I don't know what I'm doing here and so am currently treating this as a learning exercise. I definitely wouldn't trust this model to base any important decisions on but am hoping that it will be good enough to help understand what might be happening during the transition into and out of field weakening.

Also not too sure how useful it might be for MTPA as I'm not to sure how to calculate the magnitude of the reluctance torque to allow it to be included. Any ideas how to do this?

[johu]

Looking very nice.

Does this help?

grafik.png (2.82 KiB) Viewed 41239 times

Have you modeled flux linkage yet (the greek letter in the formula)?

[Pete9008]

Does this help?
grafik.png

Have you modeled flux linkage yet (the greek letter in the formula)?

Perfect, thanks!

Just realised that that is the same equation linked by Bigpie in his MTPA thread:

Still reading up at the moment. https://github.com/vedderb/bldc/ is a treasure trove, especially this commit https://github.com/vedderb/bldc/pull/179

but I hadn't twigged how to use it till your comment on flux linkage. I'm assuming that 'p' ('PP' in Bigpie's link) is the number of poles on the motor?

Afraid I wouldn't have a clue how to model the flux linkage (and even if I did it's unlikely that all the parameters needed would be available). Instead I'm just estimating it from the motor torque and power data, it's not going to be exact but hopefully good enough for now? If anyone can advise how to go about modelling it I can have a go though.

Currently the model currently just uses the

torque calc.png (1.26 KiB) Viewed 41222 times

term for torque but I'll add the reluctance term next time I work on it. It should be possible to see the effect it has by adjusting the existing Ld and Lq parameters.

Edit - just wondering how FW will affect the above equation - does it just affect the flux linkage term? If so then that implies that the reluctance torque would have a larger percentage contribution to the total torque in the FW region?

[Pete9008]

Been doing a little more reading on the maths for motors and have come to the conclusion that while the model I have at the moment is probably good enough for investigating field weakening it 's more or less useless for MTPA. Looking at a couple of posts from @bexander led to here https://se.mathworks.com/help/mcb/ref/m ... rence.html and in particular this:

At the moment my simulation just uses the first two terms of each of the first two equations (or rather it will once I add motor resistance). I don't really understand the significance of the last two terms but I'm guessing they represent cross coupling between the d and q fluxes and so model the bits needed for MTPA simulations. Without these adding the reluctance torque term is probably a waste of time.

The problem is I need to calculate Id and Iq for a given Vd and Vq so the above equations are the wrong way round. What I need is the top two equations rearranging to give Id and Iq for a given Vd and Vq. Once upon a time I might have been able to do this but currently it's beyond me so could anyone help?

Also an explanation of what the extra terms are doing would be great too, I don't like using equations without understanding them!

Edit - Looking again I am including the ωλpm term assuming that's the BEMF?

[Pete9008]

Could someone please tell me what I'm missing here?

Using this equation from above:

and the following values (combination of @bexander 's measurements and general values off the web):
p = 4
Lq=6mH
Ld=2mH
Flux linkage= 0.13Weber
Iq = 250A
Id = 0A

Gives a torque of 195Nm - seems reasonable.

With Iq or 240A and Id of -70A gives 590nM - would be nice but sounds unrealistic!

With Iq or 175A and Id of -175A gives 890nM - really???

I've seen articles suggesting a 14-40% increase in torque for MTPA but nothing like this - what am I missing?

[Pete9008]

Most recent master branch
https://github.com/jsphuebner/stm32-sin ... on-foc.cpp

You will see dependencies to other modules there, most importantly foc.cpp and picontroller.cpp which reside in libopeninv

Working on integrating this now - quick question:

In the line
uint16_t rotorAngle = Encoder::GetRotorAngle(); of PwmGeneration::CalcNextAngleSync(int dir) what are the units of rotorAngle (i.e. what is 0degree and what is 360degrees)?

Edit: and is this electrical angle or rotor angle (i.e. with or without pole pairs)?

[Pete9008]

Think I've worked out how the rotor angle code works now but in the process I came across this function:

Code: Select all

void Encoder::UpdateRotorFrequency(int callingFrequency)
{
   distance += turnsSinceLastSample;

   if ((encMode == AB) || (encMode == ABZ))
   {
      //65536 is one turn
      lastFrequency = (callingFrequency * turnsSinceLastSample) / FP_TOINT(TWO_PI);
      turnsSinceLastSample = 0;
   }
   else if ((encMode == RESOLVER) || (encMode == SPI) || (encMode == SINCOS))
   {
      int absTurns = ABS(turnsSinceLastSample);
      if (startupDelay == 0 && absTurns > STABLE_ANGLE)
      {
         lastFrequency = (callingFrequency * absTurns) / FP_TOINT(TWO_PI);
         detectedDirection = turnsSinceLastSample > 0 ? 1 : -1;
      }
      else
      {
         lastFrequency = 0;
      }
      turnsSinceLastSample = 0;
   }
}

Which is called every 10ms to calculate the rotor frequency. This is read using Encoder::GetRotorFrequency() in:

Code: Select all

void PwmGeneration::CalcNextAngleSync(int dir)
{
   if (Encoder::SeenNorthSignal())
   {
      uint16_t syncOfs = Param::GetInt(Param::syncofs);
      uint16_t rotorAngle = Encoder::GetRotorAngle();
      int syncadv = frqFiltered * Param::GetInt(Param::syncadv);
      syncadv = MAX(0, syncadv);

      //Compensate rotor movement that happened between sampling and processing
      syncOfs += FP_TOINT(dir * syncadv);

      angle = polePairRatio * rotorAngle + syncOfs;
      frq = polePairRatio * Encoder::GetRotorFrequency();

to calculate frq which is then filtered in:

Code: Select all

void PwmGeneration::Run(s32fp il1_in, s32fp il2_in, uint16_t angle_in)
{
   if (opmode == MOD_MANUAL || opmode == MOD_RUN)
   {
      static s32fp idcFiltered = 0;
      int dir = Param::GetInt(Param::dir);
      int moddedfwkp;
      int kifrqgain = Param::GetInt(Param::curkifrqgain);
      s32fp id, iq;

      Encoder::UpdateRotorAngle(dir);

      CalcNextAngleSync(dir);
      FOC::SetAngle(angle);

      frqFiltered = IIRFILTER(frqFiltered, frq, 8);

to give frqFiltered which is then used in PwmGeneration::CalcNextAngleSync(int dir) (middle code box above) to calculate the amount of syncadv to use.

If I understand the code correctly IIRFILTER(frqFiltered, frq, 8) becomes:

Code: Select all

frqFiltered = (frq + (frqFiltered*256)-frqFiltered)/256

Called every 10ms this gives a time constant of around 2.8seconds, effectively it will be using motor speed information that is seconds out of date when calculating syncadv. Not a problem in steady state but will cause odd effects when accelerating or decelerating.

I think it is quite possible that this has been causing some of the unwanted regen issues mentioned in other thread especially if the regen tends to last for 2-3 seconds.

Edit - The effect could last longer than 2-3seconds if the motor is still changing speed at the end of the period. Effectively the error in syncadv could continue until 2-3sec after the car stops decelerating/accelerating.

[johu]

Great work! I see you worked out most bits, 0°=0 digit, 360°=65536 digits

IIRFILTER is called every 113us so the time constant would be 0.113/10=0.03s

The software I published here viewtopic.php?p=43809#p43809 has the syncadv filtering disabled

[Pete9008]

Great work! I see you worked out most bits, 0°=0 digit, 360°=65536 digits

Thanks, taking a bit longer than I expected as there is a lot of code in OpenInverter and C++ is a little outside my comfort zone.

Could you confirm that the angle is the motor shaft angle not the electrical field angle (i.e. it doesn't take account of the number of poles)?

IIRFILTER is called every 113us so the time constant would be 0.113/10=0.03s

Doh! I'd got hung up on the value being updated every 10ms and forgot that the filter runs a lot faster.

Could you share the calculation you use for working out the time constant? I've always ended up doing a spreadsheet to work it out but your method looks a lot faster.

The software I published here viewtopic.php?p=43809#p43809 has the syncadv filtering disabled

OK, will stick with the code I have for now to see if I can produce the unwanted regen effect, If I can then I'll see what effect this change has

The OpenInverter code is now integrated into the simulation and is moving the motor (not properly though, just oscillating). Hopefully I just need to sort out the scaling going in and out, get all the pole pair settings so they match and get the phase offsets aligned for it to start doing something more sensible. Then it's a case of adding controls to the GUI for the important parameters before being able to run some tests.

I've had to modify the PwmGeneration source files which I was hoping to avoid doing as it makes switching to different code versions more difficult. Once it's working I'll need to go back and write a lot more test stubs to try and use the unmodified files. The main problem I had is I just can't get my head round the way you have written the anain class, the implementation is very elegant but just beyond my C++ skills. Any suggestion of how to write a test stub for it would be welcome.

[Pete9008]

Getting closer but the control loops aren't doing quite what I expect.

A question on the code, in the picontroller.cpp:

Code: Select all

int32_t PiController::Run(s32fp curVal)
{
   s32fp err = refVal - curVal;
   s32fp esumTemp = esum + err;

   int32_t y = FP_TOINT(err * kp + (esumTemp / frequency) * ki);

On the last line I believe that frequency is the call frequency, in this case around 8800Hz. If this is the case then the I term does nothing until esumTemp is greater than this. Am I missing something here?

What I've done in the past in this case is keep track of the remainder from the division and add it in on the next loop. Going to give it a try now to see if it helps.

Edit - Ignore me, I think it's fine as is!

[johu]

Could you confirm that the angle is the motor shaft angle not the electrical field angle (i.e. it doesn't take account of the number of poles)?

angle is electrical - shaft angle isn't really relevant to control.

Could you share the calculation you use for working out the time constant? I've always ended up doing a spreadsheet to work it out but your method looks a lot faster.

Hahaha, I was gonna ask you the same just took your number and divided it

[Pete9008]

Hahaha, I was gonna ask you the same just took your number and divided it

It's always bugged me having to spreadsheet it. Sure there must be a rule of thumb formulae - just don't know what it is!

The simulation seems to be working

I've ended up messing about with the OpenInverter source code more than I intended to while trying to figure out the problem so it will need a fair bit of tidying up. In then end most of it was understanding the scaling and fixed point stuff - must have tried every wrong combination before finding the right one.

First capture, just a simple transient test in manual mode with a command of 100A Q current:

This is with the PI gains from the standard parameters. After tuning it:

Next another manual mode run, 150A current but this time running till steady state is reached:

And finally a run mode run with 100% torque requested:

I think field weakening is disabled (default parameters used - need to dig into it a bit more) but it's interesting to see that the motor speed stops increasing when it gets to base frequency. Not sure how it manages to over speed slightly or why it decelerates slightly when told to coast. The later could be resistive losses in the winding, which are now included, or the fact that sync advance is running even though it's not needed (probably the next thing to have a look at).

Overall pretty happy with that. It's going to have to go on the back burner for a while now as I need to get on with the DIY that I should have been doing for the last few days rather than this but hopefully I can have a look at the field weakening in the evenings.

Edit - just noting it here so that I don't forget - the model is still missing some of the terms in the Iq and Id calculation (as described above) but is hopefully good enough for now.

[Pete9008]

Just noticed that on the second plot above the current doesn't fully return to 0 when commanded to. I'm wondering whether there is a low level error that might be fixed by carrying the remainder over in the PI calculation - need to have another look. Convinced myself it was, then convinced myself it wasn't and now back to thinking that it is

Edit - had another look and I'm now convinced it's fine as is. So that I remember when I get back to it the problem is that this section of code (in PwmGeneration::Run()) is being called which resets the I term and prevents it getting all the way back to 0 (sits at around 6A):

Code: Select all

if ((0 == frq && 0 == idref && 0 == qController.GetRef()) || initwait > 0)
      {
//         timer_disable_break_main_output(PWM_TIMER);
         dController.ResetIntegrator();
         qController.ResetIntegrator();
         fwController.ResetIntegrator();
         RunOffsetCalibration();
      }

[Pete9008]

Right, rubbish nights sleep to not safe to do DIY
On the other hand that's a good excuse to do more code (as it's just getting to the interesting bit)

Looking again at the last plot above it is obvious that:

• It's obvious that either the field weakening or MTPA algorithms are running
• My simulation isn't doing a very good job of reducing the BEMF in the presence of negative Id

So the first job today is to figure out which (wondering whether the frequency is wrong which would mess both up) and fix the model BEMF calculation.

Next is probably to disable MTPA for now, I think I've got my head round how parameters work now so will probably add one for (Ld-Lq). Then see if the model simulates FW OK.

[Bigpie]

Is your model able to run in 'run' mode and just give it a throttle position and have it calculate Id and Iq? Can we get our hands on it yet?

[Pete9008]

Is your model able to run in 'run' mode and just give it a throttle position and have it calculate Id and Iq? Can we get our hands on it yet?

Yes, you can choose which mode it runs in from the main window although I've only just started using it in run mode.

Not just yet. The model isn't at the point where I'd trust it and there are a few more controls needed in the main window before it can be useful. Also just figured out that I've definitely broken the code that updates the frequency variable that is used for FW and MPTA so that is the next job.

I've also hacked the OpenInverter code around a lot more than I wanted to in order get it working - before I release it I'd like to get it back to using unmodified OpenInverter source code (should be possible not I have a better understanding of the code - it just needs a LOT more test stubs writing ). That way you should be able to git clone whatever version of OpenInverter you want to test into the project as a subdirectory and run it.

[Pete9008]

angle is electrical - shaft angle isn't really relevant to control.

Question on the following code:

Code: Select all

#define STABLE_ANGLE      ((10 * TWO_PI) / 360)

void Encoder::UpdateRotorFrequency(int callingFrequency)
{
    distance += turnsSinceLastSample;

    int absTurns = ABS(turnsSinceLastSample);
    if (startupDelay == 0 && absTurns > STABLE_ANGLE)
    {
       lastFrequency = (callingFrequency * absTurns) / FP_TOINT(TWO_PI);
       

Is this right?
It looks like it's supposed to be only update last frequency if the angle is more than 10degrees (which makes sense) but because it's comparing to turns it needs 1820 full rotations?

Edit - please ignore, I misunderstood the code, my test harness is wrong instead.

[Pete9008]

Well that was more of a pain than it should have been!

Parameters added for Lq-Ld and flux linkage. Two problems, that took ages to find. First the fixed point representation uses a different number of decimal places in the param and foc modules meaning that paramaters can't be pulled into the foc module for the MTPA algorithm wthout conversion. Second there is a bug in the OpenInverter code MTPA algorithm. Code, with the correction, is:

Code: Select all

static const s32fp fluxLinkage = FP_FROMFLT(0.09);
static const s32fp fluxLinkage2 = FP_MUL(fluxLinkage, fluxLinkage);
//static const s32fp lqminusldSquaredBs10 = FP_FROMFLT(0.01722); //additional 10-bit left shift because otherwise it can't be represented
static const s32fp lqminusldSquaredBs10 = FP_FROMFLT(0.03444); //additional 10-bit left shift because otherwise it can't be represented
static const s32fp lqminusld = FP_FROMFLT(0.0058);

The change is the constant lqminusldSquaredBs10, the 10bit shift needs the value to be 0.03444 not 0.01722.

Screen grab for original code:

With the above fix:

And changed to use parameters gives the same result so I haven't broken anything (yet!):

This could possibly explain a few issues with reduced torque at higher revs when using MTPA!