You are here

RL measurement incorrect with factor 2/3?

32 posts / 0 new
Last post
Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17
RL measurement incorrect with factor 2/3?

Hello all,

I was measuring the resistance of a hoverboard motor using a VESC 6 MKIII. The value I get is about 140 mΩ.

I also tested the motor directly on my power supply. 0.4 V phase-phase gives 1 Amp, meaning the phase-zero resistance is about 200 mΩ (400 mV / 1 A / 2).

As the difference was quite large I started to check some things. When checking the VQ and IQ in the realtime data after the RL measurement I manually calculate close to 200 mΩ (Rs=Vq/Iq for stationary motor and current). Then in the source code I found the following in mcpwm_foc_measure_resistance:

return (voltage_avg / current_avg) * (2.0 / 3.0);

I do not think this factor 2/3 should be used here. The voltage and current in the dq domain divided by each other directly gives you the phase-zero resistance. If this factor 2/3 was not there the value would have been 210 mΩ, much closer to my measured value.

I also see the same factor used in the inductance measurement. But then, when the inductance value is used for the PMSM decoupling a factor 3/2 is applied again, making it correct there. However, as the current controller parameters use the R and L with the factor 2/3 the current controller has a control bandwidth (zero-dB crossing) of 2/3th of the setting in the calculation screen. 

I hope I am wrong and the code is correct, can somebody please explain the factor 2/3?

PS. The controller and software are really awesome, great job Benjamin!

Kind regards,

Elwin

 

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Hello Benjamin,

I checked more of the source code to get a better understanding. It looks very good, but I think there is improvement possible in some of the formulas. Here are my thoughts:

As far as I know the typical hobby BLDC motors have nearly sinusoidal back EMF. This means they can be described well using PMSM equations.

PMSM equations taken from mathworks.com/help/physmod/sps/ref/pmsm.html:

vd=Rsid+LddiddtNωiqLq,

vq=Rsiq+Lqdiqdt+Nω(idLd+ψm),

Ris the stator resistance (phase-zero). For explanation of the other symbols check the link.

 

During the resistance measurement theLqdiqdt+Nω(idLd+ψm) terms are all zero. What remains is vq=Rsiq =>Rs= v / iq

So, as far as I can see, there should be no scaling of 2/3 of the result for a PMSM.

In the code the inductance is also scaled with a factor 2/3. In the cross decoupling code actually the inductance is multiplied by 3/2, which I also did not expect (see the equations above, there is no factor 3/2 in a PMSM). However, because the inductance is reduced by a factor 2/3, the result should still be correct.

dec_vd = state_m->iq * state_m->speed_rad_s * conf_now->foc_motor_l * (3.0 / 2.0);

dec_vq = state_m->id * state_m->speed_rad_s * conf_now->foc_motor_l * (3.0 / 2.0);

 

The calculation for flux linkage also seems to try to correct for the factor 2/3, but instead of multiplying by 3/2 as in the decoupling a multiplication of 2/3 is used. This causes the estimation to incorrectly compensate for the resistive and inductive components, surely making the results invalid if the current is not small. 

*linkage = (v_mag - (2.0 / 3.0) * res * i_mag) / rad_s - (2.0 / 3.0) * i_mag * ind;

Looking back to the PMSM equations the formula used also seems a bit strange. I am not 100% sure about this but it seems more logical to me to use the following formulas instead:

vq=Rsiq+Lqdiqdt+Nω(idLd+ψm)

Lqdiqdt=0

Nω = ωe (electrical speed [rad/s])

ψm= (vq-Rsiq) / ωe -  idLd

The formula is similar in structure (except for the factor 2/3), but now the (I think) correct currents and voltages are used. I think this can improve the motor identification, especially in cases where the current is not small.

Very interested to hear your thoughts on all this. Would also be happy to help / test on my VESC.

Kind regards,

Elwin

kubark42
Offline
Last seen: 2 days 15 hours ago
VESC Free
Joined: 2020-07-17 18:46
Posts: 19

Interesting discussion. For others who might like to follow along, here are the links in the code:

https://github.com/vedderb/bldc/blob/812f5302ca7078330f4cd293abb9d33fb04...

https://github.com/vedderb/bldc/blob/3a071cee2f76db24308a9fb00b7d7482748...

-------------------

For the code authors' consideration, I believe there's also something to learn from this. The code is dependent on delicate math implementation, but little of it is documented. The upshot is that it becomes arduous to back out the theoretical underpinnings, and thus nigh impossible to evaluate and review the code. Consider this a kind coaxing for anyone implementing math to please, please, please document the crap out of it. Papers, relevant websites, the math written out in plain text, whatever it takes to make it very easy for readers to understand where the code is coming from. This will help when we're either curious to understand-- e.g. because we think we might have found an error--, or we believe we have a better algo and would like to give it a shot.

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

I think I found where the factor 2/3 comes into play. If you look at the following document about the PMSM by NXP:

https://www.nxp.com/docs/en/application-note/AN4680.pdf,

you can find that if you measure the inductance (but also resistance) like shown in figure 7 (page 9), you need to take into account a factor 2/3. This is because there is 1 inductor (or resistor) in series with 2 inductors (or resistors) in parallel.

This is however not the way the VESC should calculate R or L, because the VESC measures using vd and vq, which are not phase-phase quantities.

-------------------

Agree on the idea to better document where the used formulas come from. Thanks for the nice Github links, I did not know that was possible like that. Learn something new every day.

kubark42
Offline
Last seen: 2 days 15 hours ago
VESC Free
Joined: 2020-07-17 18:46
Posts: 19

I'm not familiar enough with PMSM motor param identification to have a well-formed opinion on the 2/3 scalar (although what you say makes a lot of sense). However I do notice a couple numerical issues in the calibration code. Both the R and the L depend on accumulators which are defined as floats. Unfortunately, the shear number of samples stresses the ability of a float to faithfully capture the summations.

In other words, a `float` can only handle a value with 2^23 = 8.4 million, which means almost 7 significant digits. That sounds like a lot, but since the R calibration uses 10,000 measurements, and the L calibration uses 20,000 measurements, then 8.4e6/10e4 = 840. So this is fewer than three significant digits for raw measurements, and barely more than two for L. I'd have to understand how much precision the ADC has to know for sure if it's an issue. However, if the `float` is indeed running out of precision, this would lead to a smaller estimate (because the `float accumulator doesn't increase but the `int` divisor  still does) which would also be consistent with your experience.

kubark42
Offline
Last seen: 2 days 15 hours ago
VESC Free
Joined: 2020-07-17 18:46
Posts: 19

I looked at the math, and I think your point rings true. If the system is measuring phase-to-phase current, then phase-zero resistance is 0.5*voltage/(phase-to-phase current). In other words, V_d/I_d. So where is the 2/3 term coming from in https://github.com/vedderb/bldc/blob/812f5302ca7078330f4cd293abb9d33fb04...?

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Yes, that is correct. However the vd vq id and iq are all related to phase values, not phase-phase. 

Phase quantities are sine waves (120 degrees apart) for a rotating motor. The (stationary) vd vq id and iq values are relate 1to1 to the peak of these sine waves. For example, if iq is 2 and id is zero (and the motor is rotating), you will have sine waves of current in the phases, with an amplitude of 2 Ampere. Similarly, a vq of 5 volt relates to sine waves (120 degrees apart) with an amplitude of 5 volt each. 

But to get back to the resistance calculation, if you look to the formulas in the Mathworks link I sent, you will find that the resistance observed in the dq domain is exactly the same as the resistance of a single phase Rs.

Note: there exists also a power invariant transform, then the relation is different. However, the vesc uses the amplitude invariant transformation, making the above explanation valid.

 

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

I did some simulations with single precision floating point numbers to check the point you raised about the precision. The way the floating point numbers are added like this introduces a small error of less then 0.1%, I would not worry about it. 

What is important is that you use sufficient current for the measurment, othewise the quantization of the adc can become a problem. But I think this is done correctly in the vesc.

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

Thanks for looking into this, it always helps to have some extra eyes on these things.

When I looked at the equations I came to the conclusion that the Clarke transform (3-phase to 2-phase) introduces the 2/3 compared to what you would measure between two motor wires, and I thought it would be useful to have the value as what you would see with a multimeter and a motor. I verified the 2/3 by running quite high current in a motor with a power supply and measure the voltage drop across it, which matched up with the 2/3 at that time. This was a long time ago though and there might have been some dead time distortion that affected my measurements, so it is possible that I got that wrong.

I also set up simulations of the whole system based on the virtual motor code contributed by axiom to make sure that everything is correct and adds up in the code. This is why the 2/3 is removed in the observer and the axis decoupling (as you saw with the axis decoupling):

https://github.com/vedderb/bldc/blob/dev_fw_5_03/mcpwm_foc.c#L3544

https://github.com/vedderb/bldc/blob/dev_fw_5_03/mcpwm_foc.c#L3725

I also verified the flux linkage calculation, and I'm fairly confident that it is correct and things add up the way they are now. That being said, all of the 2/3s in the detection, observer and decoupling code could be removed, and it would still work correctly.

One thing I'm aware of that is not correct is the time constant for the current controller parameter calculation, it will be off by a factor of 2/3. That is not too important though, it only matters that it is not too fast and not too slow and what is correct on each motor depends on the rest of the system. Most systems work well with how it is calculated now, and I have not prioritized to update it yet.

Interesting observation about the float precision, it might make sense to average over fewer samples in a few places.

 

About writing down the math, I personally often have an easier time reading code than mathematical notations (depending on both how the code and math is written obviously). Keeping two sets of the same equations updated is also a burden, and there are plenty of sources where one can read about motor math out there that should match the code more or less. Therefore I haven't spent any time making another writeup about motor math, but I encourage anyone who is interested to go through the code and document what is going in a different format. I'm sure that makes it easier to understand and follow for many people who are not like me.

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

I looked into this today with a VESC that has phase filters to work around the dead time distortion, and the 2/3 is not supposed to be there. The compensation for that was done correctly in the code at least, so everything should have been working with the exception of the PI controller time constant being calculated correctly.

I will work more on it in the weekend and report back.

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Hi Benjamin,

Thank you for the response!

Was the factor 2/3 in any way related to the dead time distortion? I do not see the link there.

--------------------------------------------------

I did the check of measuring the phase-phase resistance of my motor using my lab power supply and the result was as follows:

I = 3 A

Vph-ph = 1.61 V

Rph-ph = 1.61/3 = 0.54 Ohm

Rph-0 = Rph-ph / 2 = 0.27 Ohm.

VESC gives Rph-0 = 0.174 Ohm

However, if the factor 2/3 was not there, it would have given 0.26 Ohm, pretty much spot on compared to the 0.27 Ohm.

---------------------------------------------------

I did some more checking on the flux linkage detection and found that it is indeed correct, but depending on which option you choose.

correct: measure_linkage_foc

Because here the (incorrect) resistance is corrected with a factor 3/2:

https://github.com/vedderb/bldc/blob/c0fa922251f26e968dc64d269a6aba9878f...

(although inductance cross coupling is ignored, I think this is not a large issue)

 

not correct: measure_linkage_openloop

Here actually there is another factor 2/3 instead of 3/2 (making R actually (2/3)² too low), making the flux linkage overestimated.

https://github.com/vedderb/bldc/blob/3a071cee2f76db24308a9fb00b7d7482748...

 

Interested to hear what you find!

 

Kind regards,

Elwin

TechAUmNu
Offline
Last seen: 3 days 14 hours ago
VESC Free
Joined: 2017-09-22 01:27
Posts: 542

Can you try this with a motor with very low resistance? 

I think we also might have a problem that with low resistance motors <10mOhm. The on resistance of the mosfets will start to make quite a large impact on the result.
For the math to work should we only be wanting the motor resistance / inductance. Or do we also need to take into account the whole powerstage?
If we only want the motor resistance then maybe we need a setting in the hwconf that defines the on resistance per phase and that gets factored into the detection?

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

measure_linkage_openloop looks at both the driven and undriven flux linkage, and when it is used the undriven one is always preferred as it does not depend on the resistance and inductance. Only on motors with low inertia and high friction that does not work as the motors stops so fast that the flux linkage cannot be measured from the back emf alone, and then the driven flux linkage that has to be compensated for resistance, inductance and current is used. That being said, that is an extremely rare case and most motors measure the flux linkage in the undriven way just fine.

The problem with the formula that contains id and iq is that it assumes you know them, which you only do if you know the motor phase. All you know in openloop is the magnitude as you don't have observer parameters yet, so you cannot calculate the flux linkage properly. I know the formula is wrong, it just gets close enough in some cases.

 

About the power stage resistance, that should be included in the parameters. What matters is what the motor, power stage, connectors and cables look like together to track the motor well, including the voltage and current tolerances on the hardware. So the best resistance to use for good tracking and operation is the one you see when running the detection, not what you measure with an accurate LCR meter.

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Hi Benjamin,

Thank you for the response. Ok, that about the measure_linkage_openloop makes sense. I was not aware that the undriven value was taken.

---------------------------------------------

I was looking at the inductance measurement, which made me check out the hfi, both in code and practice. Let me start of by saying that I really like the idea and implementation, great job!

However, I think I did find a discrepancy in the math (on top of the factor 2/3), which makes the inductance measured to be lower.

In this part:

https://github.com/vedderb/bldc/blob/be439a2a55ae5b1ba8cce8ad6b0da2d8046...

                motor->m_hfi.buffer[motor->m_hfi.ind] = ((hfi_voltage / 2.0 - conf_now->foc_motor_r *
                        current_sample) / (conf_now->foc_f_sw * current_sample));

I think you try to calculate back the inductance by L = (delta_V - R*delta_i) / (delta_i / delta_t).

However, I think there are three problems here:

1. delta_V should be hfi_voltage * 2, not hfi_voltage /2 (the step in voltage is from -hfi_voltage to +hfi_voltage and visaversa, making the step 2 * hfi_voltage )

2. foc_motor_r should have the factor 3/2, unless you remove the factor 2/3 we talked about earlier. Of course this is automatically fixed if the 2/3 is removed the the R calculation as suggested.

3. I have a feeling the conf_now->foc_f_sw parameter is 1 divided by the set switching frequency. However, when measured on an external scope, one finds that the real switching frequency is actually halve of what you set. This causes a factor 2 error.

My suggestion then would be:

                motor->m_hfi.buffer[motor->m_hfi.ind] = ((hfi_voltage * 2.0 - conf_now->foc_motor_r *
                        current_sample) / (conf_now->foc_f_sw * 2 * current_sample));

Which actually doubles the measured induction.

I know my industrial servo motor specifications pretty wel, Lq = 260 uH, Ld = 220 uH. However, when measuring on the VESC it gives measurements in the range of 70 to 90 uH. Correcting for both the 2/3 factor and the above factor 2 the measurement actually would be 210 uH to 270 uH, which does make a lot of sense to me.

Really interested to hear what you think about this.

PS. Great job for making all this work!

Kind regards,

Elwin

 

kubark42
Offline
Last seen: 2 days 15 hours ago
VESC Free
Joined: 2020-07-17 18:46
Posts: 19

> Interesting observation about the float precision, it might make sense to average over fewer samples in a few places.

I made a PR to convert certain floats to doubles. The use of doubles is limited to the calibration routines, so there shouldn't be any material impact on normal performance. 

> Can you try this with a motor with very low resistance? 

I can take a look with a Geiger HPD20. It has with published specs of 3.2mOhms and 8.7uH. According to the factory, they measure resistance with a 60A current and they measure inductance at 30kHz. 

> However, if the factor 2/3 was not there, it would have given 0.26 Ohm, pretty much spot on compared to the 0.27 Ohm.

I have a similar problem with my Geiger's detected resistance. It is somewhere about 2/3 of what the manufacturer spec says it should be. This weekend I'll hook it up to my bench supply to validate.

> Correcting for both the 2/3 factor and the above factor 2 the measurement actually would be 210 uH to 270 uH, which does make a lot of sense to me.

I'm eager to see the results of this. I documented at https://vesc-project.com/node/3216 where I get a 3x smaller estimate vs the manufacturer published inductance spec.

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

I will do some more testing with the HFI and inductance measurement in the next days and report back.

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Instead of testing things with the hfi, wouldn't it make more sense to discuss the physics/math first? The behaviour can also easily be simulated. I would be happy to assist, for example in a Teams call or forum discussion, just let me know whether you are interested.

Kind regards,

Elwin

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

We can do that as I go, but I have to get into the HFI-code again and understand exactly what is going on again. There are so many details that I had to get right to make it work at all, I spent a full week implementing it doing nothing else. How do deal with sampling in V0/V7, aligning samples, making the calculation current invariant, transitioning in and out of HFI, making it run fast, using a good amount of samples, how much current to get the initial samples etc. It took me 2 hours to get a proof of concept running, but making it work well and useful was a lot of work. I have turned every stone around 10 times already.

The tricky thing with datasheets is that they don't always define resistance and inductance in the same way as the code expects it. It is quite easy to see if it is correct though with a testbench that can load a motor and has an encoder. If the inductance is wrong the observer angle will deviate from the encoder angle as the load and speed increases.

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

I agree that the details in getting everything to work right is very complicated and most of the work. However the basic concept and math is relatively easy to understand and should be the simple part to get right. I think getting these basics right should be prio 1. The good thing is you got the hard stuff all figured out already!

For example the factor 2 in the hfi injection. I think the math and physics should show whether this is needed. But, since you have gotten this far with your project, I am sure you know best on how to tackle this. 

Good luck in testing!

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Have you found time to check it out? Any comments on the formula?

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

I checked out the code a bit deeper and I think I found a flaw in my explanation.

I now see you seem to apply the hfi voltage for only halve the time between two current measurements. That would mean the voltage step is 1x hfi, not 2x hfi. Still do not see why the step would be hfi/2 though.

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

I was working on it yesterday, and that is the same conclusion I came to. All 2/3 and 3/2 are gone, and both resistance and inductance measurements are within a few percent of my LCR meter. The only confusing thing left is the division by 2 in the hfi (I didn't spend too much time in that yet though). That being said, I still have to do a lot of testing with the updated inductance measurement to make sure that motors run well under load. I will report back soon.

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

To make things clear for myself I made some hand calculations with just an inductance which shows the VESC calculation is underestimating by a factor 2. Either v_hfi should not be divided by 2, or foc_f_sw needs to be divided by 2:

 

Schermafbeelding 2021-10-14 215451.png

 

 

To also include a resistance component I used a Matlab script to simulate and come to the same conclusion:

inductance.png

So actually in the end the factor v_hfi /2 is fine, but the foc_f_sw needs to be divided by 2 (just as shown by hand calculation with just L).

Conclusions

  • The best estimator is L = (v_hfi / 2 - R * delta_i ) / ( foc_f_sw / 2 * delta_i )
  • Current VESC L estimation underestimates by a factor 2 (as I already suspected because I know my motor well, however to be absolutely sure tomorrow I will test my motor at my work on our BK Precision 891 LCR meter)
  • When time constant L/R becomes very low, the inductance measurement is off. Too much transient behaviour during cycle: these motors therefore need higher switching frequency anyhow.
  • I do not see how your LCR measurement can match VESC measurement with the current formula.

 

 

Matlab code for anyone interested:

TAU = 1 ./ logspace(0 , 4 , 20);
for jmodel = 1:length(TAU) %s time constant of electrical system (tau = L/R)
    
    tau = TAU( jmodel );
    L = 100e-6;  %[H]   Inductance (in the direction the hfi is injecting)
    R = L / tau; %[Ohm] Resistance (in the direction the hfi is injecting)
    
    v_hfi = 1; %[V] voltage level for hfi
    v_bus = 48; %[V]
    
    foc_f_sw = 30e3; %[Hz]
    
    t = linspace( 0 , 2/foc_f_sw , 1e5); %[s] Time vector
    
    v = zeros(size(t));
%     v(t< 1/foc_f_sw) =  v_hfi; %[V] input voltage vector. v_hfi for half the time of foc_f_sw.
    v(t< (1/foc_f_sw) * v_hfi / v_bus) =  v_bus; %[V] input voltage vector. Takes into account bus voltage.
    
    sys = tf( 1 , [L R]); %Transer function of system I/V = 1/(Ls+R)
    
    i = lsim( sys , v , t ); %[A] Simulate current
    
%     figure(1)
%     subplot( 2, 1 , 1)
%     plot( t*1e6 , v , 'linewidth' , 2 )
%     hold all
%     ylabel('Voltage [V]')
%     xlabel('Time [\mus]')
%     grid on
%     subplot( 2, 1 , 2)
%     plot( t*1e6 , i , 'linewidth' , 2 )
%     hold all
%     xlabel('Time [μs]')
%     ylabel('Current [A]')
%     grid on
    
    di = i(end) - i(1);

    L_est1(jmodel)      = (v_hfi   - R * di) / (foc_f_sw   * di );
    L_est2(jmodel)      = (v_hfi/2 - R * di) / (foc_f_sw/2 * di );
    L_est_vesc(jmodel)  = (v_hfi/2 - R * di) / (foc_f_sw   * di );
    
end

figure
semilogx( TAU([1 end]) , [L L] *1e6 , 'linewidth' , 2 , 'Color' , 'Black')
hold all
semilogx( TAU , L_est1 *1e6 , 'linewidth' , 2 )
semilogx( TAU , L_est2 *1e6, 'linewidth' , 2 )
semilogx( TAU , L_est_vesc *1e6 , 'linewidth' , 2 )
grid on
xlabel('Electrical time contant (L/R) [s]')
legend( '"Real" L value: 100 μH' , 'Option1: L = (v_hfi - R * di) / (foc_f_sw * di )',  'Proposal: L = (v_hfi/2 - R * di) / (foc_f_sw/2 * di )' , 'VESC now: (v_hfi/2 - R * di) / (foc_f_sw * di )' )
ylabel('Inductance [μH]')
title('hfi indutance estimation simulation (foc\_f\_sw = 30 kHz)')

Jfriesen
Offline
Last seen: 2 days 18 hours ago
VESC Free
Joined: 2017-10-24 22:49
Posts: 6

Nice work here! Ben and I have been investigating this and I did some more work with your matlab script. Thanks! 

First step was to get the positive and negative pulses being applied:

photo_2021-10-15_12-44-25.jpg

Looking here we can see the real issue is the exponential rolloff since the samples are being taken at the center of the off time here:

photo_2021-10-15_12-48-22.jpg

so I made a very ugly drawing:

photo_2021-10-15_17-01-02.jpg

And tried to find a closed form solution using the solution to the differential equation of an RL circuit as described here

Unfortunately that was an absolute mess and I am fairly certain no easy solution exists. So instead I opted for a fixed point solution. We can start with our existing estimate for R and L and solve for ix1 and ix2 in the picture above. Then I solved the response between ix1 and ix2 with the applied negative bus voltage for L. Iterate through this a few times updating L and hopefully we converge. I got lazy and used symbolic to solve the exponential equation it and the result is this: 

 L = (R*t_duty)/log((-V_bus - Ix2*R)/(-V_bus - Ix1*R))

Then the total fixed point thing looks like this: 

L_guess = v_hfi / (foc_f_sw * di );
i_high = mean(i(end/10:2*end/10:end));
i_low = mean(i(1:2*end/10:end));
for xxx = 1:10        
   Ix1 = i_high * exp(-R/L_guess*t_off/2);
   Ix2 = i_low * exp(R/L_guess*t_off/2);
   L_guess = -(R*t_on)/log((-v_bus - Ix2*R)/(-v_bus - Ix1*R))
end
L_est2(jmodel)      = -(R*t_on)/log((-v_bus - Ix2*R)/(-v_bus - Ix1*R));

The mean part is just picking off the points and obviously isn't needed here in this world of perfect precision. 

The results at high and low duty: 

photo_2021-10-15_16-59-23.jpg

photo_2021-10-15_17-00-46.jpg

Can see the compensated line is falling dead on hooray! 

Full code below: 

%pkg load control
close all
TAU = 1 ./ logspace(0.1 , 4 , 20);
for jmodel = 1:length(TAU) %s time constant of electrical system (tau = L/R)
    
    tau = TAU( jmodel );
    L = 100e-6;  %[H]   Inductance (in the direction the hfi is injecting)
    R = L / tau; %[Ohm] Resistance (in the direction the hfi is injecting)
    
    v_hfi = 3; %[V] voltage level for hfi (no impact on result in linear model)
    v_bus = 48; %[V]
    
    foc_f_sw = 10e3; %[Hz]
    
    t = linspace( 0 , 500/foc_f_sw , 4e5); %[s] Time vector
    
    v = zeros(size(t));
    for iii = 1:250
    v((t > 0.5 * (1/foc_f_sw*(4*iii - 3 - v_hfi / v_bus))) & ( t < 0.5 * (1/foc_f_sw*(4*iii - 3 + v_hfi / v_bus)))) =  v_bus; %[V] input voltage vector. v_hfi for half the time of foc_f_sw.
    v((t > 0.5 * (1/foc_f_sw*(4*iii - 1 - v_hfi / v_bus))) & ( t < 0.5 * (1/foc_f_sw*(4*iii - 1 + v_hfi / v_bus)))) =  -v_bus; %[V] input voltage vector. v_hfi for half the time of foc_f_sw.
    end    
    
    sys = ss(tf( 1 , [L R])); %Transer function of system I/V = 1/(Ls+R)
    x0 = 0;
    for jj = 1:10
        [~, ~, x] = lsim( sys , v , t, x0 ); %[A] Simulate current
        x0 = x(end);
    end
    
    t = linspace( 0 , 10/foc_f_sw , 4e5); %[s] Time vector
    
    v = zeros(size(t));
    for iii = 1:5
    v((t > 0.5 * (1/foc_f_sw*(4*iii - 3 - v_hfi / v_bus))) & ( t < 0.5 * (1/foc_f_sw*(4*iii - 3 + v_hfi / v_bus)))) =  v_bus; %[V] input voltage vector. v_hfi for half the time of foc_f_sw.
    v((t > 0.5 * (1/foc_f_sw*(4*iii - 1 - v_hfi / v_bus))) & ( t < 0.5 * (1/foc_f_sw*(4*iii - 1 + v_hfi / v_bus)))) =  -v_bus; %[V] input voltage vector. v_hfi for half the time of foc_f_sw.
    end    
    
    sys = ss(tf( 1 , [L R])); %Transer function of system I/V = 1/(Ls+R)
    [i, tout, x] = lsim( sys , v , t, x0 ); %[A] Simulate current
    
     figure(1)
     subplot( 2, 1 , 1)
     plot( t*1e6 , v , 'linewidth' , 2 )
     hold all
     ylabel('Voltage [V]')
     xlabel('Time [\mus]')
     grid on
     subplot( 2, 1 , 2)
     plot( t*1e6 , i , 'linewidth' , 2 )
     hold all
     xlabel('Time [μs]')
     ylabel('Current [A]')
     grid on
     
    i(end/10:2*end/10:end)
    di = mean(i(end/10:2*end/10:end)) - mean(i(1:2*end/10:end));
    avgi = (mean(i(end/10:2*end/10:end)) - mean(i(1:2*end/10:end)))/2;

    L_est1(jmodel)      = (v_hfi - 0 * R * avgi) / (foc_f_sw * di );
    t_on = 1/foc_f_sw*v_hfi/v_bus;
    t_off = 1/foc_f_sw*(1-v_hfi/v_bus);
    
    L_guess = v_hfi / (foc_f_sw * di );
    i_high = mean(i(end/10:2*end/10:end));
    i_low = mean(i(1:2*end/10:end));
    for xxx = 1:10        
        Ix1 = i_high * exp(-R/L_guess*t_off/2);
        Ix2 = i_low * exp(R/L_guess*t_off/2);
        L_guess = -(R*t_on)/log((-v_bus - Ix2*R)/(-v_bus - Ix1*R))
    end
    L_est2(jmodel)      = -(R*t_on)/log((-v_bus - Ix2*R)/(-v_bus - Ix1*R));  
    
end

figure
semilogx( TAU([1 end]) , [L L] *1e6 , 'linewidth' , 2 , 'Color' , 'Black')
set(gca, "fontsize", 14)
hold all
semilogx( TAU , L_est1 *1e6, 'linewidth' , 2 )
semilogx( TAU , L_est2 *1e6, 'linewidth' , 2 )

grid on
xlabel('Electrical time contant (L/R) [s]')

legend('"Real" L value: 100 μH','v__hfi  / (foc_f_sw * di )','compensated')

ylabel('Inductance [μH]')
title('hfi indutance estimation simulation (foc\_f\_sw = 30 kHz)')

 

Jeff

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Wow, you did some nice work there! Seems to be a very good estimation like that. However, we should not forget that even if you can estimate L correctly, having these enormous current ripples is pretty bad for your performance and efficiency. Motors which such a low time constant really beg for much higher switching frequencies. That being said, having a good estimation for L seems pretty helpful to me.

 

As promised, today I finished fully measuring another of my industrial motors on a BK Precision 891 LCR meter and the VESC 6 MK III. The results are as follows:

 

BK Precision 891 LCR meter:

Rphase-zero = 258 mOhm (Raw measurement 516 mOhm phase-phase)

Ld = 167 μH (Raw measurement between phase 1 and 2,3 together 250 μH with rotor aligned to D axis. Ld is 2/3 of this measurement, see NXP PMSM document)

Lq = 212 μH (Raw measurement between phase 1 and 2,3 together 318 μH with rotor aligned to Q axis. Lq is 2/3 of this measurement, see NXP PMSM document)

Lavg = (Lq + Ld)/2 = 189.5 μH

LdLqdiff = 45 μH

 

VESC:

Rphase-zero = 165,66 mΩ (= 248.5 if the 2/3 is removed)

Lavg = 66,14 µH  (= 99.21 µH if the 2/3 is removed, 198.4 µH if the L estimation is also fixed)

ld_lq_diff: 7.33 uH  (= 11 µH if the 2/3 is removed, 22 µH if the L estimation is also fixed)

Note how everything will be very close if the corrections are made, expect for the ld_lq difference (unless you define it as offset from the average to both sides, but that seems to not be the variable name). I Think you need to double the ld_lq_diff measurement and do the other fixes described earlier and everything seems to be very close.

I had a question about the HFI rotor tracking. Does it assume the higher inductance is the D axis or the Q axis? I see for my motor the Q axis has a higher induction, but this can also be the other way around, see for example the rotor configurations in: https://en.engineering-solutions.ru/motorcontrol/pmsm/.

Kind regards,

Elwin

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

In my previous post I already mentioned the ld_lq difference should be increased by a factor of 2. I checked the code and found out why this happens:

The FFT of the signal is taken and then divided by the number of samples. Then correctly the abs of this value is taken. However, what seems to be forgotten is the FFT of a frequency component gives 2 components in the frequency domain, one at the positive frequency, and one at the negative frequency. For any real input signal (without any imaginary part) the component at the negative frequency is the same value as the positive frequency. To get back the amplitude of the time domain signal one has to add the negative and positive frequency components. To achieve this one can simply multiply all (except for the DC component) results by 2. This factor 2 seems to be missing in the code, explaining why the ld_lq difference is a factor 2 too low. 

Afterwards there correctly is a factor 2 to go from zero-peak to peak-peak.

Example figures of what happens in case the amplitude is 1:

 

cos.pngfft.png

Matlab code which shows the effect:

x = 1:32;
y = cos(2 * x/32 *2*pi);

figure
plot( x , y)
xlabel('Measurement point')
ylabel('Amplitude')
title('cos(2 * x/32 *2*pi)')

n = length(y);
Y = fft(y) / n;

Y2 = fftshift(Y);
fshift = (-n/2:n/2-1); 

figure
plot( fshift , abs(Y2) )
xlabel('Bin')
ylabel('Magnitude')
xlim( [-10 10])
title('Frequency domain')

 

 

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

Everything looks correct on my end now, the only mystery left is the /2 in the HFI code. For now I will just ignore it and continue testing.

This is a known inductor Jeff made that we used for testing. It shows the same value (within a few percent) in both my LCR meter and the VESC detection.

20211015_113456.jpg

I also tested one of my outrunners and carefully measured the lowest and highest inductance with the LCR meter to get the average inductance and ld_lq_diff (at 1 kHz, 10 kHz). The VESC is very close to this one too (half of what is measured between 2 phases, which should correspond to one winding).

After that I did some testing on a known difficult motor: the T-motor U15 with a large propeller

20211015_174606.jpg

It worked better than ever now at up 250 phase amps (iq). Before I had to tweak the inductance, but now I just took the default values and increased the observer gain by a factor of 4.

I got some promising results with another difficult motor too, but on this one I had to halve the observer gain.

Next I will connect an encoder to my dual motor rig and see if the observer angle drifts away from the encoder angle as I increase the load and speed - that should give an indication if the inductance is correct.

Everything is pushed on github and the latest VESC Tool beta version with the updated firmwares is updated here: https://vesc-project.com/node/2859

Please give it a try and let me know how it works!

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Hi Benjamin,

Great job, looks very promising! Good to hear that things are working even better, and with less tweaking!

 

I have a few small remaining remarks on the code changes:

-----------------------------------------------------

You removed the scaling here,

https://github.com/vedderb/bldc/blob/d8a99fd0f1d065eb75ecdb8b742d1883e46...

but now your real voltage is a factor 1.5 higher than what you use in the calculation. The scaling should be like you have in other parts of the code:

https://github.com/vedderb/bldc/blob/d8a99fd0f1d065eb75ecdb8b742d1883e46...

(although I really don't understand the 2/3 that is used, but I checked and the outcome with the SVM function gives the correct amplitudes for the phases)

I think this factor 1.5 is why the results are now very close, without fixing the L calculation function.

-----------------------------------------------------

https://github.com/vedderb/bldc/blob/d8a99fd0f1d065eb75ecdb8b742d1883e46...

Power loss in alpha-beta and dq domain are given by 1.5*i*i*R, not 2*i*i*R.  

For a given i, the 3 phase currents are sine waves with i as amplitude. the RMS value of these waves is equal to i / sqrt(2). The power loss for each of these is (i / sqrt(2))^2* R. Since you have 3 phases you end up with Ploss = 1.5*i*i*R.

-----------------------------------------------------

These lines need scaling, or your alpha-beta voltages can be higher than you can make without overmodulation (you correcty removed the 2/3, but you need scaling there):

https://github.com/vedderb/bldc/blob/d8a99fd0f1d065eb75ecdb8b742d1883e46...

https://github.com/vedderb/bldc/blob/d8a99fd0f1d065eb75ecdb8b742d1883e46...

The correct scaling between duty cycle and dq or alpha-beta voltages is:  duty * mc_interface_get_input_voltage_filtered()  * ONE_BY_SQRT3.

This is because the largest sine waves you can make have phase-phase voltage of Vbus. However, the dq/alpha beta are related to phase-zero quantities, which are 1/sqrt(3) the amplitude of the phase-phase quantity.

-----------------------------------------------------

Kind regards,

Elwin

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

I did the math, and you are correct about the loss calculation. It did not add up with the measurements I made, and it turns out that the power error was due to dead-time distortion in the bus current estimation. This has been an issue for a long time, but now it is improved a bit.

The HFI max voltage should be fixed now as well. I agree that the 2/3 should be there and the division by 2 should be removed, but I simply don't get correct measurements if I do that on any of the hardware versions. I tested the inductor and all my motors with the detection and LCR meter, and this is what gives the most correct results by far under a range of input voltages, switching frequencies and HFI voltages. For now I will just ignore that and move on. At least the other things in the code are consistent now with what you would read in a datasheet of a motor, and the inductance value works much better. I will still do some testing with an encoder on my dual motor rig and confirm that it tracks well in the next days, hopefully that reveals something else.

The test builds are currently building, I will update the link in the other post in a few minutes.

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

I agree that it is best to move on. Happy to see accurate measurement results. 

benjamin
Offline
Last seen: 11 hours 32 min ago
Joined: 2016-12-26 15:20
Posts: 422

I installed an encoder on my dual motor rig, and it turns out that it works better with the (slightly higher) induction measured by correct formulas:

code.png

When the inductance is measured like this, the phase error by the observer when applying a load becomes smaller compared to the same test with the other inductance. That indicates that it is closer to the truth, even if it is a bit higher than what I measured with the LCR meter (although the other formula gave lower inductance than the LCR meter).

I will continue testing in the next days and report back!

Elwin
Offline
Last seen: 1 hour 37 min ago
VESC Free
Joined: 2021-09-30 16:41
Posts: 17

Sounds good. Good luck with your testing!