import numpy as np
from SciDataTool.Functions.Plot.plot_3D import plot_3D
from ...Classes.Electrical import Electrical
from ...Classes.EEC_PMSM import EEC_PMSM
from ...Classes.Simu1 import Simu1
from ...Classes.OPdq import OPdq
from ...Classes.InputCurrent import InputCurrent
[docs]def comp_MTPA(
machine,
LUT,
N0_min,
N0_max,
Nspeed,
Ntorque,
I_max,
U_max,
Rs=None,
Tsta=20,
n_Id=100,
n_Iq=100,
is_plot=False,
):
"""Return the Air Gap Surface Force
Parameters
----------
machine : machine
A machine object
LUT: LUTdq
A look up table object
Tem_vect: ndarray
Array of requested torque level
N0_vect: ndarray
Array of requested speeds
I_max: float
Maximum current [Arms]
U_max: float
Maximum voltage [Vrms]
Rs: float
Stator phase resistance to enforce [Ohm]
Tsta: float
Stator widning temperature to enforce [deg]
n_Id: int
Desciretization number to interpolate along d-axis
n_Iq: int
Desciretization number to interpolate along q-axis
Returns
-------
OP_matrix_MTPA: ndarray
OP_matrix resulting from MTPA strategy
U_MTPA: ndarray
Output voltage [Vrms]
I_MTPA: ndarray
Output current [Arms]
"""
# Speed vector
N0_vect = np.linspace(N0_min, N0_max, Nspeed)
# Maximum load vector
if Ntorque == 1:
Tem_vect = np.array([])
else:
Tem_vect = np.linspace(0, 1, Ntorque)
# Pole pair number
p = machine.get_pole_pair_number()
# Stator winding number of phases
qs = machine.stator.winding.qs
OP_matrix = LUT.get_OP_array("N0", "Id", "Iq")
# Get Id_min, Id_max, Iq_min, Iq_max from OP_matrix
Id_min = OP_matrix[:, 1].min()
Id_max = OP_matrix[:, 1].max()
Iq_min = OP_matrix[:, 2].min()
Iq_max = OP_matrix[:, 2].max()
Id, Iq = np.meshgrid(
np.linspace(Id_min, Id_max, n_Id), np.linspace(Iq_min, Iq_max, n_Iq)
)
Id, Iq = Id.ravel(), Iq.ravel()
Imax_interp = np.sqrt(Id ** 2 + Iq ** 2)
simu_MTPA = Simu1(name=machine.name + "_comp_MTPA", machine=machine)
# Definition of the input
OP_ref = OPdq(N0=OP_matrix[0, 0], Id_ref=OP_matrix[0, 1], Iq_ref=OP_matrix[0, 2])
simu_MTPA.input = InputCurrent(OP=OP_ref)
# Stator winding resistance [Ohm]
ecc = EEC_PMSM(R1=Rs) # If Rs None => Compute
simu_MTPA.elec = Electrical(LUT_enforced=LUT, eec=ecc, Tsta=Tsta)
# Interpolate stator winding flux in dqh frame for all Id/Iq
eec_param = simu_MTPA.elec.eec.comp_parameters(
machine=machine, OP=OP_ref, Tsta=simu_MTPA.elec.Tsta, Id_array=Id, Iq_array=Iq
)
# Compute torque
Tem_sync, Tem_rel = simu_MTPA.elec.eec.comp_torque_sync_rel(eec_param, qs, p)
Tem_interp = Tem_sync + Tem_rel
# Init OP_matrix
OP_matrix_MTPA = np.zeros((Nspeed, Ntorque, 4))
U_MTPA = np.zeros((Nspeed, Ntorque, 3))
I_MTPA = np.zeros((Nspeed, Ntorque, 3))
for ii, N0 in enumerate(N0_vect):
print("Speed " + str(ii + 1) + "/" + str(Nspeed))
# Update operating point
OP_ref.N0 = N0
OP_ref.felec = None
eec_param = simu_MTPA.elec.eec.comp_parameters(
machine,
OP=OP_ref,
Tsta=simu_MTPA.elec.Tsta,
Trot=simu_MTPA.elec.Trot,
Id_array=Id,
Iq_array=Iq,
)
# Calculate voltage
out_dict = simu_MTPA.elec.eec.solve(eec_param)
U_max_interp = np.sqrt(out_dict["Ud"] ** 2 + out_dict["Uq"] ** 2)
# Finding indices of operating points satisfying Vmax and I_max voltage and torque limitations
j0 = np.logical_and(U_max_interp <= U_max, Imax_interp <= np.abs(I_max))
# Finding index of operating point giving maximum positive torque among feasible operating points
jmax = np.argmax(Tem_interp[j0])
# Maximum torque achieved for N0
Tem_max = Tem_interp[j0][jmax]
# Store values in MTPA
OP_matrix_MTPA[ii, -1, 0] = N0
OP_matrix_MTPA[ii, -1, 1] = Id[j0][jmax]
OP_matrix_MTPA[ii, -1, 2] = Iq[j0][jmax]
OP_matrix_MTPA[ii, -1, 3] = Tem_max
U_MTPA[ii, -1, 0] = out_dict["Ud"][j0][jmax]
U_MTPA[ii, -1, 1] = out_dict["Uq"][j0][jmax]
U_MTPA[ii, -1, 2] = U_max_interp[j0][jmax]
I_MTPA[ii, -1, 0] = OP_matrix_MTPA[ii, -1, 1]
I_MTPA[ii, -1, 1] = OP_matrix_MTPA[ii, -1, 2]
I_MTPA[ii, -1, 2] = Imax_interp[j0][jmax]
for kk, Tem_rate in enumerate(Tem_vect[:-1]):
if Tem_rate == 0:
Tem_k = 0.01 * Tem_max
# Finding indices of operating points satisfying Vmax voltage for no torque production
j0 = np.logical_and(U_max_interp <= U_max, np.abs(Tem_interp) < Tem_k)
# Finding index of operating point for lowest current
jmax = np.argmin(np.abs(Imax_interp[j0]))
else:
# Finding indices of operating points satisfying maximum voltage and torque level
j0 = np.logical_and(
U_max_interp <= U_max, Tem_interp >= Tem_rate * Tem_max
)
# Finding index of operating points with lowest current among feasible operating points
jmax = np.argmin(Imax_interp[j0])
# Store values in MTPA
OP_matrix_MTPA[ii, kk, 0] = N0
OP_matrix_MTPA[ii, kk, 1] = Id[j0][jmax]
OP_matrix_MTPA[ii, kk, 2] = Iq[j0][jmax]
OP_matrix_MTPA[ii, kk, 3] = Tem_interp[j0][jmax]
U_MTPA[ii, kk, 0] = out_dict["Ud"][j0][jmax]
U_MTPA[ii, kk, 1] = out_dict["Uq"][j0][jmax]
U_MTPA[ii, kk, 2] = U_max_interp[j0][jmax]
I_MTPA[ii, kk, 0] = OP_matrix_MTPA[ii, kk, 1]
I_MTPA[ii, kk, 1] = OP_matrix_MTPA[ii, kk, 2]
I_MTPA[ii, kk, 2] = Imax_interp[j0][jmax]
if is_plot:
# Init plot map
dict_map = {
"Xdata": Id.reshape((n_Iq, n_Id))[0, :],
"Ydata": Iq.reshape((n_Iq, n_Id))[:, 0],
"xlabel": "d-axis current [Arms]",
"ylabel": "q-axis current [Arms]",
"type_plot": "pcolor",
"is_contour": True,
}
# Plot torque maps
plot_3D(
Zdata=Tem_interp.reshape((n_Iq, n_Id)).T,
zlabel="Average Torque [N.m]",
title="Torque map in dq plane",
**dict_map,
)
# plot_3D(
# Zdata=Id.reshape((n_Iq, n_Id)).T,
# zlabel="Average Torque [N.m]",
# title="Torque map in dq plane",
# # save_path=join(save_path, name + "_torque_map.png"),
# **dict_map,
# )
# plot_3D(
# Zdata=Iq.reshape((n_Iq, n_Id)).T,
# zlabel="Average Torque [N.m]",
# title="Torque map in dq plane",
# # save_path=join(save_path, name + "_torque_map.png"),
# **dict_map,
# )
# plt.contour(
# dict_map["Xdata"],
# dict_map["Ydata"],
# U_max_interp.reshape((n_Iq, n_Id)),
# colors="blue",
# linewidths=0.8,
# )
plot_3D(
Zdata=Tem_sync.reshape((n_Iq, n_Id)).T,
zlabel="Synchrnous Torque [N.m]",
title="Torque map in dq plane",
**dict_map,
)
plot_3D(
Zdata=Tem_rel.reshape((n_Iq, n_Id)).T,
zlabel="Reluctant Torque [N.m]",
title="Torque map in dq plane",
**dict_map,
)
# Plot Phi_d map
plot_3D(
Zdata=eec_param["Phid"].reshape((n_Iq, n_Id)).T,
zlabel="$\Phi_d$ [Wb]",
title="Flux linkage map in dq plane (d-axis)",
**dict_map,
)
# Plot Phi_q map
plot_3D(
Zdata=eec_param["Phiq"].reshape((n_Iq, n_Id)).T,
zlabel="$\Phi_q$ [Wb]",
title="Flux linkage map in dq plane (q-axis)",
**dict_map,
)
return OP_matrix_MTPA, U_MTPA, I_MTPA
