IPM Motor Simulation and FOC Software

[Pete9008]

And one more interesting plot for now. It's the same as the last one but with MTPA disabled by setting Ld and Lq the same so all Iq:

Bit early to draw any conclusion but it would appear that it gives much better acceleration. It would aslo appear that I need to sort out the sync advance code as it is causing an error between the currents seen by the motor and controller. The currents as seen by the controller are in the Debug window and look correct. The motor currents are seen in the Motor Data window and it can be seen that there is an angle error increasing with speed which is causing a negative Id. This is effectively acting as FW and allowing the motor to exceed base speed!

[Pete9008]

And in case it's usefull to anyone - to add MTPA parameters the following changes are needed to the code (only the changed sections of the files included):

param_prj.h

Code: Select all

//Next param id (increase when adding new parameter!): 141
//Next value Id: 2049
/*              category     name         unit       min     max     default id */

#define MOTOR_PARAMETERS_COMMON \
    PARAM_ENTRY(CAT_MOTOR,   polepairs,   "",        1,      16,     2,      32  ) \
    PARAM_ENTRY(CAT_MOTOR,   respolepairs,"",        1,      16,     1,      93  ) \
    PARAM_ENTRY(CAT_MOTOR,   sincosofs,   "dig",     1,      4096,   2048,   131 ) \
    PARAM_ENTRY(CAT_MOTOR,   encmode,     ENCMODES,  0,      5,      0,      75  ) \
    PARAM_ENTRY(CAT_MOTOR,   fmax,        "Hz",      21,     1000,   200,    9   ) \
    PARAM_ENTRY(CAT_MOTOR,   numimp,      "ppr",     8,      8192,   60,     15  ) \
    PARAM_ENTRY(CAT_MOTOR,   dirchrpm,    "rpm",     0,      20000,  100,    87  ) \
    PARAM_ENTRY(CAT_MOTOR,   dirmode,     DIRMODES,  0,      4,      1,      95  ) \
    PARAM_ENTRY(CAT_MOTOR,   snsm,        SNS_M,     12,     23,     12,     46  )

#define MOTOR_PARAMETERS_SINE \
    PARAM_ENTRY(CAT_MOTOR,   boost,       "dig",     0,      37813,  1700,   1   ) \
    PARAM_ENTRY(CAT_MOTOR,   fweak,       "Hz",      0,      1000,   90,     2   ) \
    PARAM_ENTRY(CAT_MOTOR,   fweakstrt,   "Hz",      0,      1000,   400,    134 ) \
    PARAM_ENTRY(CAT_MOTOR,   fconst,      "Hz",      0,      1000,   180,    99  ) \
    PARAM_ENTRY(CAT_MOTOR,   udcnom,      "V",       0,      1000,   0,      78  ) \
    PARAM_ENTRY(CAT_MOTOR,   fslipmin,    "Hz",      0.3,    10,     1,      37  ) \
    PARAM_ENTRY(CAT_MOTOR,   fslipmax,    "Hz",      0.3,    10,     3,      33  ) \
    PARAM_ENTRY(CAT_MOTOR,   fslipconstmax,"Hz",     0,      10,     5,      100 )

#define MOTOR_PARAMETERS_FOC \
    PARAM_ENTRY(CAT_MOTOR,   curkp,       "",        0,      20000,  32,     107 ) \
    PARAM_ENTRY(CAT_MOTOR,   curki,       "",        0,      100000, 20000,  108 ) \
    PARAM_ENTRY(CAT_MOTOR,   curkifrqgain,"dig/Hz",  0,      1000,   50,     120 ) \
    PARAM_ENTRY(CAT_MOTOR,   fwkp,        "",        -10000, 0,      -100,   118 ) \
    PARAM_ENTRY(CAT_MOTOR,   ffwstart,    "Hz",      0,      1000,   200,    136 ) \
    PARAM_ENTRY(CAT_MOTOR,   syncofs,     "dig",     0,      65535,  0,      70  ) \
    PARAM_ENTRY(CAT_MOTOR,   syncadv,     "dig/Hz",  0,      65535,  10,     133 ) \
    PARAM_ENTRY(CAT_MOTOR,   lqminusld,   "mH",      0,      1000,   5.8,    139 ) \
    PARAM_ENTRY(CAT_MOTOR,   fluxlinkage, "mWeber",  0,      1000,   90,     140 )

foc.cpp - comment out the lines shown and add the params.h include (at the top of the file)

Code: Select all

#define CST_DIGITS 15
#include "my_fp.h"
#include "my_math.h"
#include "foc.h"
#include "sine_core.h"
#include "params.h"

#define SQRT3 FP_FROMFLT(1.732050807568877293527446315059)
#define R1 FP_FROMFLT(0.03)
#define S1 FP_FROMFLT(0.15)
#define R2 FP_FROMFLT(0.5)
#define S2 FP_FROMFLT(0.5)
#define S3 FP_FROMFLT(1)
#define RADSTART(x) x < R1 ? S1 : (x < R2 ? S2 : S3)

//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);
static const u32fp sqrt3 = SQRT3;

and change the Mtpa function to:

Code: Select all

void FOC::Mtpa(int32_t is, int32_t& idref, int32_t& iqref)
{
   s32fp fluxLinkage = (Param::Get(Param::fluxlinkage)<<10)/1000;//additional 10-bit left shift due to CST_DIGITS difference between here and param module
   s32fp fluxLinkage2 = FP_MUL(fluxLinkage, fluxLinkage);
   s32fp lqminusld = (Param::Get(Param::lqminusld)<<10)/1000;
   s32fp lqminusldSquaredBs10 = FP_MUL(lqminusld<<10, lqminusld);//additional 10-bit left shift because otherwise it can't be represented
   int32_t isSquared = is * is;
   int32_t sign = is < 0 ? -1 : 1;
   if(lqminusld != 0)
   {
       s32fp term1 = fpsqrt(fluxLinkage2 + ((lqminusldSquaredBs10 * isSquared) >> 10));
       idref = FP_TOINT(FP_DIV(fluxLinkage - term1, lqminusld));
   }
   else
       idref = 0;
   iqref = sign * (int32_t)sqrt(isSquared - idref * idref);
}

Edit - Error in param_prj.h file above fixed, range for lqminusld and fluxlinkage should be 0 -1000

[Bigpie]

I did a build with and went for a drive. Results was pretty underivable

Code: Select all

static const s32fp lqminusldSquaredBs10 = FP_FROMFLT(0.03444); //additional 10-bit left shift because otherwise it can't be represented

I've attached logs incase it's of interest. logs (6) is with the FP_FROMFLT(0.03444); logs (7) is with as it was.

[Pete9008]

That was quick!

is "underivable " good or bad?

Will have a look at the logs now.

[Pete9008]

Looking at the data I'm guessing bad

Interesting, on the simulation it is definitely better and in terms of the maths it makes sense (unless I'm missing something?). It is possible that other parameters have been tuned with this error there to minimise its effect?

[Pete9008]

Been playing around with the effect of the MTPA algorithm and either there is something wrong with my simulation or something wrong with the algorithm. It seems to add much higher levels of Id than I would expect.

Going to spreadsheet the equations later to see if they make more sense that way but for now be very careful doing anything based on the above!

[Bigpie]

Bad very little torque. The flux linkage and d q difference is almost certainly not correct for my motor though

[johu]

Hmm, interesting! I fixed the MTPA formula at some point, so your change plus a change in the formula itself and got less performance and strange behaviour on cruise control (or probably rather at low load). So I reverted the change and never looked at it again so far.

Let me know if I can help getting your code on github

[Pete9008]

Hmm, interesting! I fixed the MTPA formula at some point, so your change plus a change in the formula itself and got less performance and strange behaviour on cruise control (or probably rather at low load). So I reverted the change and never looked at it again so far.

Let me know if I can help getting your code on github

I'd not looked at the formulae but just have and it's not quite what I was expecting based on the mathworks reference equations:

Are these the equation that should be implemented or are there other versions?

Was it the lqminusldSquaredBs10 value change that you tried before, what was the formulae change?

Re git, would it be possible to post a zip of the qt project here on the forum (once I've got it so it will work with unmodified OpenInverter builds)?

[Pete9008]

I think the MTPA algorithm in the code is:

Code: Select all

Id = (1/(Lq-Ld)) x (λ - sqrt(λ2 + (Is2 x (Lq-Ld)2))))

While the Mathworks one rearranges to:

Code: Select all

Id = (1/(4 x (Lq-Ld))) x (λ - sqrt(λ2 + (8 x Is2 x (Lq-Ld)2))))

So similar but not the same, any thoughts?

[Pete9008]

Some more plots (let me know if this is using too much forum storage space!) :

First 4 captures with the current MTPA algorithm (top equation in last post).
First is with MTPA Off

Second is with the Controller Lq-Ld matched to the Motor Lq-Ld

Third is with the Controller Lq-Ld lower than the Motor Lq-Ld

Fourth is with the Controller Lq-Ld higher than the Motor Lq-Ld

Next 4 captures with the revised MTPA algorithm (lower equation in last post).
First is with MTPA Off

Second is with the Controller Lq-Ld matched to the Motor Lq-Ld

Third is with the Controller Lq-Ld lower than the Motor Lq-Ld

Fourth is with the Controller Lq-Ld higher than the Motor Lq-Ld

Comments welcome! If I had to chose one it is the new matched, especially compared to the old matched, higher acceleration to a faster speed (note that this is all without field weakening). The old matched is particularly bad but it's interesting that the old plot with a high value of Lq-Ld is almost as good as the new matched one.

If trying this in a car (with the parameter control added) I'd start with a guessed flux linkage (note that this is in mWeber in the new param_prj.h file) somewhere in the right ballpark. Set Lq-Ld to zero (to disable MTPA) and carefully try it (based on the simulation torque should be low). Then increase Lq-Ld by 1mH at a time and see whether things improve. Based on forum posts I would guess something around 4-6mH might be best. Note that an incorrect value of Lq-Ld will have an effect on field weakening too so best stay well away from that region while testing! Please note that I'm not suggesting trying it yet (I wouldn't) but just in case anyone feels the need...

Edit - For those that noticed there are now more plots - I realised that the old low Lq-Ld was actually with MTPA turned off (Lq-Ld = 0) so have rerun the low plots with Lq-Ld = 1. Old low plots now included as MTPA Off.

And for completeness here is the revised Mtpa function (old code left for reference but commented out):

Code: Select all

void FOC::Mtpa(int32_t is, int32_t& idref, int32_t& iqref)
{
   s32fp fluxLinkage = (Param::Get(Param::fluxlinkage)<<10)/1000;//additional 10-bit left shift due to CST_DIGITS difference between here and param module
   s32fp fluxLinkage2 = FP_MUL(fluxLinkage, fluxLinkage);
   s32fp lqminusld = (Param::Get(Param::lqminusld)<<10)/1000;
   s32fp lqminusldSquaredBs10 = FP_MUL(lqminusld<<10, lqminusld);//additional 10-bit left shift because otherwise it can't be represented
   int32_t isSquared = is * is;
   int32_t sign = is < 0 ? -1 : 1;
   if(lqminusld != 0)
   {
       s32fp term1 = fpsqrt(fluxLinkage2 + ((lqminusldSquaredBs10 * isSquared) >> 7));//Note - only 7bit shift to give x8
       idref = FP_TOINT(FP_DIV(fluxLinkage - term1, (lqminusld<<2)));
//       s32fp term1 = fpsqrt(fluxLinkage2 + ((lqminusldSquaredBs10 * isSquared) >> 10));//Note - only 7bit shift to give x8
//       idref = FP_TOINT(FP_DIV(fluxLinkage - term1, (lqminusld)));
   }
   else
       idref = 0;
   iqref = sign * (int32_t)sqrt(isSquared - idref * idref);
}

[johu]

I like the direction in which this is going! Also this adds the possibility to only run q-current for non-IPM motors
I think ultimately I'd add a function like FOC::SetMotorConstants() that pre-computes various terms

[Ev8]

Great work guys

[johu]

Are these the equation that should be implemented or are there other versions?

Was it the lqminusldSquaredBs10 value change that you tried before, what was the formulae change?

I think I'd found a more simplified version from TI buried somewhere in the FOC thread I suppose. I guess they all do the same though.
Look here: https://github.com/jsphuebner/libopenin ... 888b8a3ebf

Re git, would it be possible to post a zip of the qt project here on the forum (once I've got it so it will work with unmodified OpenInverter builds)?

yes of course

[Pete9008]

I think ultimately I'd add a function like FOC::SetMotorConstants() that pre-computes various terms

Good idea, that would be the best way to do it. It's neater and avoids including the param module into foc which presents plenty of scope for confusion on decimal places. It would also allow the conversion to be done using FRAC_DIGITS and CST_DIGITS (both of which would then be visible) rather than hard coded (which would break if either of the _DIGITS defines were changes. I won't do it now because it would break compatibility of the simulator with old code builds.

Look here: https://github.com/jsphuebner/libopenin ... 888b8a3ebf

Thanks, that gives me confidence I'm moving in the right direction. It's almost the same as the new code but missing the 1/4 scaling factor which means Id would be 4x too big. That's might explain why Id seems so big in a lot of captures.

More plots! The last thing to sort out before I can look at the regen issue is to implement SyncAdv support. I've modded my model so that it can give out an offset position (controlled by the SyncDelay contol on the main window in us). The idea is that this can be tuned to compensate for the SyncAdv in the OpenInverter code and then allow the effect of discrepancies to be investigated.

First plot is with SyncAdv set to 0 (to disable it) and then SyncDel tuned till the controller (Debug window) and motor (Motor window) views of Iq and Id match :

It's interesting that 113us is needed to get the match. This is exactly one control loop so I'm guessing that the model isn't quite doing the loop value update at the right moment. Not going to worry about it for now though.

Next plot is with the OpenInverter SyncAdv increased to the standard value of 10 and the model left at the zero position (actually 113 till I fix it):

You can clearly see the effect this has on Iq and Id as seen by the motor - a gentle shift in value with increasing speed.

And finally with SyncDel tuned to get back to correct operation:

Whats really interesting here is that it needs 280us, after allowing for the 113us error this is a advance of 167us. I was expecting a figure nearer 50us. All I can think is the the OpenInverter code is somehow giving the same 113us loop delay as the model but can't think why.

Any ideas? Is 167us possible or is the standard value of 10 for SyncAdv a bit too high?

I think this is good enough for now though so the plan is to put the old MTPA code back in and go looking for the unexpected regen!

[Bigpie]

I'll try a build with the scaled MTPA algorithm, I've been trying to measure my motor parameters, all I've got so far the the bemf constant 23.72 vrmf (span it with a drill and got 59.61Hz with a 60V peak to peak) Not managed to measure d or q inductance with any confidence yet though.

[Pete9008]

I'll try a build with the scaled MTPA algorithm, I've been trying to measure my motor parameters, all I've got so far the the bemf constant 23.72 vrmf (span it with a drill and got 59.61Hz with a 60V peak to peak) Not managed to measure d or q inductance with any confidence yet though.

As long as you are somewhere close with the bemf constant it should be OK. Since all you need for MTPA is Lq-Ld I'd just start at 1mH and increase in mH steps. From others posts I'd guess it would be around 4mH?? Make sure you include all the MTPA code mods though otherwise results could be very unpredictable.

I may have modelled the unwanted regen

Have a look at the following plots and see what you think:

First plot, field weakening turned off - limits at base frequency:

FW enabled, ramps up to high speed (higher than I would expect given fmax is set to 10krpm??):

Accelerating to a speed below base and the reducing torque request to zero - coasts:

And finally accelerate to above base freq and then come off the throttle - heavy deceleration!:

Does this look like the problem people have been experiencing?

[Ev8]

Yes that looks like it!

[Pete9008]

Couple more plots.

First with SyncAdv disabled:

Still have the same issues, strong regen.

Second with the new MTPA algorithm:

More interesting, still get the regen, although not as violent, but also now have issues with the FW (does not spin up as far) and less powerfull acceleration (previous plots were using 70% throttle, this used 100% and still slower.

So questions:
1) Has the FW been tuned to work best with the old MPTA or is there a problem with the new MTPA?
2) Do we believe the acceleration of the old MPTA? It uses very strong Id and I'm not sure I believe the torque calculation from the mathworks equations in that case )see earlier post).
3) If the regen isn't due to SyncAdv or due to the old MPTA bugs then what is causing it?

[Pete9008]

And one more plot. Just noticed the other graph window (one I haven't been showing as it's not normally that interesting), the one that shows the voltages that the controller is generating:

Whats, very interesting is that there is one blip at the point of deceleration, zooming in:

So there is a momentary blip (as the control loops stabilise) but no overall change to the commanded Vq and Vd. Don't think that this is what I expected, would expect a change.

Going to have a look at the frequency filtering mentioned in another thread to see what effect reducing the time constant on that has.

[Bigpie]

Done a build with the new MTPA, not had chance to go testing yet. Hopefully get out tomorrow. I've not experienced the regen issue personally, but I don't use lift off regen only when brake pressed. So that's not really something I can test.

[Pete9008]

Reducing the filter order from 8 to 2 gave no change so it's not that.

I don't think that this is an algorithm problem, I think it's a control loop problem. If you look at the above graphs the controller is trying to drive at max voltage on all three phases. This point also coincides with where the FW stops working properly. Essentially it has run out of control parameters, all it can change is the phase but that is only 1 parameter. It's trying to control 2, Iq and Id. Controlling 2 outputs with 1 input isn't possible. Instead of limiting the PWM outputs we need to start reducing the torque request to maintain some headroom on the pwm amplitude in order to maintain control of both amplitude and phase.

Does that make sense or am I just spouting nonsense? (just thinking out load really)

[Pete9008]

Done a build with the new MTPA, not had chance to go testing yet. Hopefully get out tomorrow. I've not experienced the regen issue personally, but I don't use lift off regen only when brake pressed. So that's not really something I can test.

No problem, thanks for trying it out. Does this have the new parameters for flux linkage and Lq-Ld?

[Pete9008]

I don't think that this is an algorithm problem, I think it's a control loop problem. If you look at the above graphs the controller is trying to drive at max voltage on all three phases. This point also coincides with where the FW stops working properly. Essentially it has run out of control parameters, all it can change is the phase but that is only 1 parameter. It's trying to control 2, Iq and Id. Controlling 2 outputs with 1 input isn't possible. Instead of limiting the PWM outputs we need to start reducing the torque request to maintain some headroom on the pwm amplitude in order to maintain control of both amplitude and phase.

Fairly sure this is what is going on. The current firmware doesn't start adding FW current until Iq fails to meet its target value. The trouble is that this condition doesn't happen until after the voltage/pwm is already maxed out and at this point we have lost one of our control variables. Instead I think the FW controller needs to have a target of keeping the PWM output at say 97% with a controller output that can only increase -Id current (i.e. if the pwm output is low it can't force it to increase).

@johu , what do you think?

[Pete9008]

Just tried the suggested control loop (in a very hacky way) and it looks promising:

Loop stability needs a bit of work but it does keep the voltage low enough to maintain full control.

Edit - This is with MRPA disabled.

Thoughts