Extension of Electrical Equivalent Circuits (EEC) to time harmonics

[jlebesnerais]

The EEC implemented in Pyleecan (voltage to current transfer function) deal with fundamental values of voltage / currents, but when supplied with PWM or when introducing back emf harmonics the equivalent circuit does not work. Extending EEC to time harmonics requires first to introduce skin effect and inductance factors for R, L elements depending on shape of conductor (round, square), without proximity effects for the moment. Skin effects would also allow to improve winding loss calculations.

[BonneelP]

Hello all,

After further discussions with my colleagues here are our new proposal for the Electric module organization:

• We want to introduce a new class "ELUT" for Electric Look-Up Table (new property of Electrical module). The point is to store all the meaningful output of the Electric module (EEC parameters, inductance tables...) and the method to interpolate them on the proper frequencies. It would contains:
  • For PMSM Phid(Id, Iq), Phiq(Id, Iq)
  • For Asynchronous machine: Lm(Im),
  • For both resistances, leakage inductance
• For the EEC, it would still have the comp_parameter method but it would check if the parameter can be interpolated first from the ELUT, else it would compute it and store it in the ELUT. Then the module can still be used without an existing ELUT, the point is to store these outputs in a convenient way for further simulations.
• The Electrical class would still handle the workflow of the model. In particular, it would compute/solve the EEC for each frequency.

Best regards, Pierre

[Sijie-NI]

Hello,

I want to calculate magnetic fluxes Phid, Phiq (over time) from three-phase frame (PhiA, PhiB, PhiC), with Park transformation. Then, I load these fluxes in SciDataTool variable form (DataTime). Could I output directly the phase, module, frequency of their harmonics by "get_harmonics"?

Best regards Sijie

def comp_EMF_harmonics(Phi_A, Phi_B, Phi_C,f, fs):

Phi_d=sqrt(2/3)*(Phi_A*cos(delta)+Phi_B*cos(delta-2*pi/3)+Phi_C*cos(delta+2*pi/3))
Phi_q=sqrt(2/3)*(-Phi_A*sin(delta)-Phi_B*sin(delta-2*pi/3)-Phi_C*sin(delta+2*pi/3))
Phi_h=sqrt(2/3)*1/2*(Phi_A+Phi_B+Phi_C)

time_axis = DataLinspace(
    name="time",
    unit="s",
    initial=0,
    final=(len(Phi_d) - 1) / fs,
    number=len(Phi_d),
    include_endpoint=True,
)
Phi_d_data = DataTime(
    name="Phi_d",
    symbol="Phi_d",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_d,
)

Phi_q_data = DataTime(
    name="Phi_q",
    symbol="Phi_q",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_q,
)
Phi_h_data = DataTime(
    name="Phi_h",
    symbol="Phi_h",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_h,
)

# Number of harmonics
N_harm=2000

Phi_d_fft=Phi_d_data.get_harmonics(N_harm)
Phi_q_fft=Phi_q_data.get_harmonics(N_harm)
Phi_h_fft=Phi_h_data.get_harmonics(N_harm)

v_d=-2*pi*f*Phi_q_fft.values*1j
v_q=2*pi*f*Phi_d_fft.values*1j

return v_d, v_q

[Sijie-NI]

Hello,

New release about comp_BEMF_harmonics:

Python code

def comp_BEMF_harmonics(Phi_A,Phi_B,Phi_C, delta, time): """ Compute the back electromotive force harmonics from magnet fluxes (PMSM)

Parameters
----------
Phi_A: Magnetic flux of phase A
Phi_B: Magnetic flux of phase B
Phi_C: Magnetic flux of phase C
delta: Rotor angular position
time: time vector

"""

#Park transformation (keep the amplitude factor=2/3)
Phi_d=2/3*(Phi_A*cos(delta)+Phi_B*cos(delta-2*pi/3)+Phi_C*cos(delta+2*pi/3))
Phi_q=2/3*(-Phi_A*sin(delta)-Phi_B*sin(delta-2*pi/3)-Phi_C*sin(delta+2*pi/3))
Phi_h=2/3*1/2*(Phi_A+Phi_B+Phi_C)

# Create time vector in form of DataLinspace
time_axis = DataLinspace(
    name="time",
    unit="s",
    initial=0,
    final=time[-1],
    number=len(time),
    include_endpoint=True,
)

# Load Phi into DataTime
Phi_A_data = DataTime(
    name="Phi_A",
    symbol="Phi_A",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_A,
)

Phi_B_data = DataTime(
    name="Phi_B",
    symbol="Phi_B",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_B,
)
Phi_C_data = DataTime(
    name="Phi_C",
    symbol="Phi_C",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_C,
)

Phi_d_data = DataTime(
    name="Phi_d",
    symbol="Phi_d",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_d,
)

Phi_q_data = DataTime(
    name="Phi_q",
    symbol="Phi_q",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_q,
)
Phi_h_data = DataTime(
    name="Phi_h",
    symbol="Phi_h",
    unit="Wb",
    normalizations=None,
    axes=[time_axis],
    values=Phi_h,
)

# Phi_q_data.plot_2D_Data("time")
# Phi_q_data.plot_2D_Data("freqs", type_plot='curve')
# plt.show()

# Calculate FFT for Phi on the dq0 frame
d = Phi_d_data.get_along("freqs")
freqs_d = d['freqs']
complx_d = d['Phi_d']

q = Phi_q_data.get_along("freqs")
freqs_q = q['freqs']
complx_q = q['Phi_q']

h = Phi_h_data.get_along("freqs")
freqs_h = h['freqs']
complx_h = h['Phi_h']

# Calculate back-emf (E) on dq0 frame
E_d=-2*pi*freqs_q*complx_q+2*pi*freqs_d*complx_d*1j
E_q=2*pi*freqs_d*complx_d+2*pi*freqs_q*complx_q*1j
E_h=2*pi*freqs_h*complx_h*1j

return E_d, E_q, E_h, freqs_d, freqs_q, freqs_h

Does it correspond with your needs? especially for the input variables of the function? Best regards, Sijie

[SebGue]

Hello Pierre,

I wasn't aware of the ELUT thing until now. Have you been working on this issue recently? I use something similar, i.e. a "MagLUT" and some extended EEC with table interpolation. So we may share some ideas if you like.

BR S

[jlebesnerais]

Heb Sebastian, Yes we are working on Electrical Look Up Table right now. The idea is to include magnetizing inductances function of currents (Lm(I) for SCIM, Phid and Phiq(Id,Iq) for PMSM) and have an Equivalent Electrical Circuit to feed the models with voltage instead of currents. The ELUT makes the V to I transfer function while the MagLUT makes the I to B transfer function. What have you done so far on your side ?

[SebGue]

Hello,

is there already a branch on your work? I got a MagLUT that is essentially a special XOutput of a mulit simulation over a grid of some reasonable d- and q-axis current densities. (I use current densities to be independent of the actual winding number of turns.) So this is my 'fluxlinkage over current' table. I also use this data to calculate iron losses, etc. The interpolated table data are then used in some eec to calculate the machine operational parameters, e.g. to get the well known efficiency maps. So at the moment I only calculate effective values. There is are no harmonics incluced. On the other hand I have no explicit transfer function V to I. I use some optimization to solve the eec. So I can impose/constrain V and/or I.

Since you may not only want to solve for effective values but also harmonics, how do you plan to realize the transfert function? Will it be time domain or frequency domain? For the magnetizing functions, they should be very similar to my MagLUT.

If you got something ready, I'm very interested to see :-)

BR S

[EmileDvs]

Hello @SebGue,

We are currently working on the EEC branch of EOMYS-public: https://github.com/EOMYS-Public/pyleecan/tree/EEC This branch also introduces many improvements on how to handle time/angle axes + periodicities, because we are also working on EEC for on-load SCIM (inclduing non-zero slip).

Concerning what we are doing for EEC of PMSM, the first step is to implement a method in OutMag to calculate phi_d and phi_q from phi_wind for a given Id/Iq point and return it as a DataFreq.

Then we'll build the ELUT the same as you did, i.e. a variable load simulation for given Id/Iq list and in the end we apply the method to calculate Phi_d/Phi_q from all phi_wind in the output_list.

Finally we get Ld/Lq function of id/iq by derivation of Phi_d and Phi_q, and back-emf for Id=Iq=0 Arms.

From this we can run the EEC and deduce current fundamental and harmonics, e.g. due to PWM harmonics.

Do you plan to share your code on how to interpolate in between a grid of Id/Iq currents and how to solve EEC using optimization ?

I ask specifically about the optimization because we are also working on using the ELUT to calculate maximum torque per ampere (MTPA) strategy, and I could find a solution but it would be cleaner by solving an optimization problem and I couldn't manage to make it work.

We are going to make a clean pull request of our branch and explain all changes.

Best regards, Emile

[SebGue]

Hello @EmileDvs,

This is a really huge topic with lots of aspects. In general I think it would make sense to not having 2 approaches for the same 'issue'. So I'm totally fine if we bring our ideas/code together. The question is how we could do that efficiently.

For the EEC of PMSM, why DataFreq? Is it for the rotor angle?

Is @Superomeg4 also involved in this. I started some discussion with him on something related, i.e. some idea how to hopefully improve the interpolation on the ELUT. Nevertheless I can share the grid interpolation and also the optimization. Just tell me where/how. Do you have any praticular reason using Ld/Lq instead of the flux linkage directly? (Btw. for completeness you got to use the coupling inductance Ldq = Lqd too.) In my opinion the equation system will be more simple using flux linkage. Otherwise you got to differentiate between absolute and differential inductances also.

One feature of my MagLUT, I missed to tell but that seems important to me, is that it is independent of some machine properties, e.g. actual length, number of turns, ... So I can reuse the table if I change lenght or winding without costly recalculation of the table.

For the optimization, maybe you can already give me your EEC/MTPA code and I adapt my solution to it?

I'm really interessted in your ideas / implementation.

Best regards, Sebastian

[EmileDvs]

Hello @SebGue,

Thanks for your answer, we are indeed several persons working on this topic, @Superomeg4 included.

For me the most efficient way to merge our approaches is:

1. on our side we continue to work on EEC branch (cf link on previous post), which contain all the changes that we have already done/planned to do. As said before, we will summarize these changes in a global issue (#452, to be completed) for more clarity.
2. we create a branch on Eomys-Public based on EEC branch, let's say EEC_merge for example
3. you implement the new features you want to share in EEC_branch and if it requires some rework/merge with the existing models in EEC branch we'll work on this together. Maybe the best here is to present the main changes you want to add in a new issue
4. make a pull request of EEC_merge to pyleecan/master that link the global issue #452 and the new issue you created with your improvements

Please let us know if you are ok with this,

Best regards, Emile

[EmileDvs]

@SebGue, here are also some answers to your other questions/remarks:

For the EEC of PMSM, why DataFreq? Is it for the rotor angle?

Yes I see your point, to be honest I am still not sure for now if we should save in time or frequency domain in ELUT. It's mainly a question of optimization, i.e. in which domain we will need more these quantities. And since all EEC for now are in frequency domain, it seems more optimal to me to store data as DataFreq. Anyway there is an equivalence between DataTime and DataFreq. Besides, we have introduced new normalizations in SciDataTool that enable to switch any quantity from time to rotor angle vectors, accounting for speed (non necessarily constant) and for initial rotor position.

Do you have any praticular reason using Ld/Lq instead of the flux linkage directly? (Btw. for completeness you got to use the coupling inductance Ldq = Lqd too.) In my opinion the equation system will be more simple using flux linkage. Otherwise you got to differentiate between absolute and differential inductances also.

Yes the aim is to calculate PWM currents from PWM voltages at variable frequency so the inductance is more suitable to get variable frequency impedance. You are completely right about calculating coupling inductances and the assumption Ldq=Lqd looks fine to me. We will probably use the dq model developed equations 4/5/6:

( reference https://www.researchgate.net/publication/335828937_Structure-borne_Noise_at_PWM_Excitation_using_an_Extended_Field_Reconstruction_Method_and_Modal_Decomposition)

One feature of my MagLUT, I missed to tell but that seems important to me, is that it is independent of some machine properties, e.g. actual length, number of turns, ... So I can reuse the table if I change lenght or winding without costly recalculation of the table.

Really interesting it is definitely a pyleecan feature that would be powerful for optimization/parameter sweep.

For the optimization, maybe you can already give me your EEC/MTPA code and I adapt my solution to it?

The EEC code is already in EEC branch, however MTPA is a script that I built before the work we did on the branch so I have to adapt it the new changes. I'll mention you whenever it is ready :)

Best regards, Emile

[SebGue]

Hello @EmileDvs ,

I think we will get wrong results, if we use this equations system. The reason is, that they use the differential (or incremental) inductance in the 'w_e * L' terms, while it should be the absolute inductance there. Do you have any particular reason to use this system? Otherwise we should derive the equations ourself and to our needs.

Best regards, Sebastian

[SebGue]

@EmileDvs I'm okay with your approach. But I got one wish :-) Maybe you can do a mini roadmap of who is doing which part and when. This is since I don't know where to start. E.g. the EEC seems to be the same as (long) before (except for skin eff.).

[EmileDvs]

Hello @EmileDvs ,

I think we will get wrong results, if we use this equations system. The reason is, that they use the differential (or incremental) inductance in the 'w_e * L' terms, while it should be the absolute inductance there. Do you have any particular reason to use this system? Otherwise we should derive the equations ourself and to our needs.

Best regards, Sebastian

My mistake, these equations using incremental inductance are for transient calculation of id/iq, this is not the proper reference. We have already developed the equations somewhere, let us check and we'll come back to you.

[EmileDvs]

@EmileDvs I'm okay with your approach. But I got one wish :-) Maybe you can do a mini roadmap of who is doing which part and when. This is since I don't know where to start. E.g. the EEC seems to be the same as (long) before (except for skin eff.).

Yes the idea of #452 is to explain what changes we are doing on our side. For the when, we are actively developing PWM currents calculation right now. You are right concerning PMSM EEC, the changes are not committed at the moment so there is no changes compared with long before.

Most of EEC changes are for SCIM, however we didn't describe it in a proper issue, so I am writing it.

The other changes (big structural changes) are listed in #452 .

As I suggested before, you can create a new issue and a new branch on EOMYS-public in which you add the models you have already developed and that you want to share. The issue will be added to #452.

[EmileDvs]

Concerning the calculation of currents harmonics due to PWM, we finally chose in PR #464 to implement a simplified analytical model (EEC_ANL object) that neglects phase resistance, such as:

1.id_pwm = vd_pwm/(j*Ld*w)

2. iq_pwm = vq_pwm/(j*Lq*w)

with w=electrical pulsation of the considered harmonic of pwm, with w >> w_elec (the fundamental electrical pulsation) to ensure that resistance is lower than inductance impedance. It is basically the same as the voltage equations below with Rs=0, ws=0 and dphi/dt= j*L*w*I:

Picture originally posted by @Sijie-NI in https://github.com/Eomys/pyleecan/issues/413#issuecomment-893239904

The already existing EEC_PMSM does not include yet voltage and current harmonics, it is still work in progress and will be developed in a future PR.

[EmileDvs]

issue solved with PR #464 and PR #481