RL measurement incorrect with factor 2/3?

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.

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.

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...?

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.

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.

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?

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.

> 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.

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

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.

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.

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:

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:

so I made a very ugly drawing:

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 i_x1 and i_x2 in the picture above. Then I solved the response between i_x1 and i_x2 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: 

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)')

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.

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

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!

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.

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:

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!

Well spotted. Curious to see how much of a difference it makes to the accuracy of the observer.

@Benjamin, in a few hours I'm going to try tuning a HPD20, which is a 20kW continuous 42-pole motor, optimized for peak power at ~2500rpm. It should be similar to the T-Motor U15, only maybe even more cantankerous. I'd love to try these code fixes, if they're pushed somewhere where I could give them a shot. Are they available online yet? I checked the `dev_fw_5_03` branch but I didn't see them there.

You can find the beta firmware on the first message on this topic : https://www.vesc-project.com/node/2859

Using the 2021-10-16 beta firmware, on my HPD20 the autodetection gives 3.50mH @30kHz for a single phase, which is not too far off from the published 8.7uH @30kHz for phase-to-phase. +-20% seems really decent for this measurement, and might be all we can hope for out of what is otherwise a motor controller and not a sensitive LCR meter.

The motor seems to start more easily, but that could also be because of other changes or even because I'm using PPM instead of keyboard control.

Only up to 180A so far, but running strong with the auto-detected settings.

That's good news, I hope you will be able to drive this motor at full power 

Did you have to tweak the observer gain ?

As near as I can tell[*], this is rock solid! I didn't have to do any tuning at all for the system to start and run impeccibaly. I got to 300A[**] with nary an instability in the motor current. Now my next concern is overheating, as even with the ESC in the propwash-- which is about as good a fan as one can get-- I was only able to do 180A[***] sustained.

[*] I do have an issue where the battery bus current estimation is off by 20%. I have validated the true current with a very accurate current shut, as well as the power supply's built-in ammeter.
[**] Only 250A once correcting for the VESC's incorrect current estimation.
[***] 150A once correcting for the VESC's incorrect current estimation.

For the overheating concern, you can try to attach the ESC to an aluminum plate, as large as possible, to maximize the heat exchange surface and put some good thermal paste between

@Elwin, I did not measure the phase currents. I can probably jury-rig something together to shunt across some of my 1.0Ohm power resistors. The motor probably won't run right, but it doesn't need to for a 20% error validation.

It would be ideal if I could set my power supply to some low voltage-- e.g. 12V-- and turn the FET completely on, running the current through the 1.0Ohm resistor. That would give ideal results.

@NexusG, in the above case I had attached some finned heat sinks with thermal paste. It was just thrown together so it wasn't optimal, but it does somewhat bound the problem.

Increasing mass will only buy time until thermal overload. That's not always a terrible thing-- in my case I've only got 3 minutes of batteries at full throttle-- however, in general it's not a steady-state solution and in flight applications extra mass is to be avoided.

I won't say I'm a pro in heat sinks, but I did my early graduate work in heat transfer so I've got a foundation to build on. Experience tells me it's going to be hard to triple-- or even quadruple-- heat dissipation without getting creative. Better would be to reduce the amount of heat which needs dissipating.

I have pushed the latest changes and updated the build here:

https://vesc-project.com/node/2859

After more testing with encoders it seems that the tracking under load is best with the new inductance measurement mode.

I also updated the observer to account for motor saliency when the Ld-Lq setting is non-zero, similar to how it was done in a paper titled A Nonlinear Observer for High-Performance Sensorless Speed Control of IPMSM Drive. That seems to work quite well.

Regarding the input current measurement of the 100/250, it is much closer than 20% for me with the latest firmware that takes the dead-time distortion into account (I got around 5% error). That was with a 40 uH motor though, so maybe the current ripple of a 4 uH motor will cause more errors in the measurement.

Small observation during testing: Running the measure_ind with very high duty cycle gives strange results for me:

I gave it a try on the VESC6 MK5, and 0.89 works for me. What was the input voltage and switching frequency?

I think this is quite unwise for some of the changed we/you made: not everything follows directly from the math. Most clear example for me is the factor 2 you changed to 4 here:

If someone is up for writing down the math down I fully support that. For now I have put a link in the source code to your post at that place to make it less confusing. I'm just saying that I haven't given that priority myself and that it isn't that high up on my priority list.

I probably won't be using VESC much more anyhow so I cannot be bothered too much.

I hope the reason for this isn't the response you got here. I really appreciate your input and assistance, it was super helpful. If there is anything I can do to encourage you to change your mind please let me know.

That makes sense, sounds like you are going in a slightly different direction. That being said, I want to move the VESC code into that direction too (the goal is to be good at all motor control applications). The main reason I have focused on robust sensorless control so far is that most people who use the VESC benefit the most from that in their applications. I would be happy to collaborate on high performance features if you are interested, and add support for all of that in VESC Tool as we go. One step in that direction would be to make some hardware with EtherCAT for high speed low latency communication - the STM32F4 should work fine for that. To make the control loop run faster there are a few things that can be done too, like optimizing the code and having a mode that skips all the observer and sensor calculations.

Some final remarks about the inductance measurement:

1) Having the resistance in the calculation made almost no difference for most motors, and Jeffs simulation shows that it can make things worse in some situations. It would be good to implement his fixed point solution if there are motors that don't work too well with the current implementation.

2) The main reason, as far as I can tell, that the inductance is over-estimated is that the current measurement amplifiers and hardware introduce some latency. When the current is sampled it has not risen to the full value yet, resulting in a slightly lower di that results in a larger inductance estimate. The reason that this works better with the observer is that the current will have the same latency there, making the apparent inductance higher. As far as I can tell the latency and the apparent higher inductance should cancel out with the observer, but maybe that isn't true in all operating points. For calculating kp and ki the actual inductance might be a better starting point, but it should not matter that much.

If someone is up for writing down the math down I fully support that

@benjamin, I have opened a GitHub issue as a proposed avenue toward better documenting the code:

https://github.com/vedderb/bldc/issues/367

I know at first glance they might not seem related, but it will be easier to propose the code comments if there is a public-facing searchable deeply-linked anchor to those changes. In particular, this indirect context helps document the kind of inscrutable `2 *` --> `4 *` changes that @elwin mentioned above. Maybe this could be done a the same time as https://github.com/vedderb/bldc/issues/362?

@Elwin, I'm a controls engineer and I would just love to know more about the features in your controller. As you say, VESC is the more mature platform and so I'm very keen to see it grow in my applications. (You've already made a huge input, because with the sys ident fixes running my huge motor just went from daunting to trivial.) In particular, I'm a data hound and am interested in gleaning as much as possible from the data feed. Having full resolution into the controller's sensor data means that it becomes straightforward to debug any odd behavior, as we can determinstically disambiguate the code from the sensors. I can see tremendous benefit to VESC to have a similar feature.

Correct, with 1.0 the svm function gives the full length of one of the six vectors it can produce. The sqrt(3) / 2 comes from the vectors in between. It is quite easy to see if you draw a line in the middle between two of the vectors:

The 2/3 then comes from the Clarke transform:

https://electronics.stackexchange.com/questions/98614/why-there-is-a-fac...

Defining it this way means that 100% modulation outputs the full input voltage in the abc-frame, which then introduces the 2/3 when moving to the alpha-beta-frame. It can be done in different ways of course, but I think this is as good as any other. Let me know what you think.