ac_motor.py

import numpy as np
from math import pi, sqrt
from scipy.integrate import odeint, complex_ode
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

filenames = [
    "dc_test_blue_motor_2019-04-05.json",
#    "blocked_rotor_test_blue_motor_2019-04-05.json"
]

class ACMotor():
    """
    Models an induction motor based on Eq 10 in [1].
    [1] https://pdfs.semanticscholar.org/4770/15e472da4c2e05e9ff8c1b921c76a938f786.pdf
    """

    def __init__(self,
                 stator_resistance, # aka r_s, [Ohm]
                 stator_inductance, # aka l_s, [Henry]
                 rotor_resistance, # aka r_r [Ohm]
                 rotor_inductance, # aka l_r [Henry]
                 magnetization_inductance # aka l_m [Henry]:
                ):
        self.magnetization_inductance = magnetization_inductance
        self.leakage_factor = 1 - magnetization_inductance**2 / (rotor_inductance * stator_inductance) # aka sigma [unitless]
        self.coupling_factor = magnetization_inductance / rotor_inductance              # aka k_r [unitless]
        self.r_sigma = stator_resistance + self.coupling_factor**2 * rotor_resistance   # [Ohm]
        self.tau_rotor = rotor_inductance / rotor_resistance                            # [s]
        self.tau_stator_prime = self.leakage_factor * stator_inductance / self.r_sigma  # [s]

        # Assigned in run():
        self.stator_voltage = None
        self.omega_stator = None
        self.omega_rotor = None

    def system_function(self, t, y):
        rotor_flux, stator_current = y # aka Phi_r, i_s [Wb, A], complex numbers

        # [1] Eq 10a
        dstator_current_dt = ((
                -1j * self.omega_stator * self.tau_stator_prime * stator_current
                - self.coupling_factor / (self.r_sigma * self.tau_rotor) * (1j*self.omega_rotor * self.tau_rotor - 1) * rotor_flux
                + 1 / self.r_sigma * self.stator_voltage
            ) - stator_current) / self.tau_stator_prime

        # [1] Eq 10b
        drotor_flux_dt = ((
                -1j * (self.omega_stator - self.omega_rotor) * self.tau_rotor * rotor_flux
                + self.magnetization_inductance * stator_current
            ) - rotor_flux) / self.tau_rotor

        return [drotor_flux_dt, dstator_current_dt]

    def run(self, time_series, voltage, omega_stator, omega_rotor):
        self.stator_voltage = voltage
        self.omega_stator = omega_stator
        self.omega_rotor = omega_rotor

        ys = np.nan * np.ones([2, time_series.size], dtype=np.complex)
        ys[:, 0] = [0.0, 0.0] # initial condition

        r = complex_ode(self.system_function)
        r.set_initial_value(ys[:, 0], time_series[0])

        for i, t in enumerate(time_series[1:]):
            if r.successful():
                ys[:, i + 1] = r.integrate(t)

        return ys

# load test data from disk

realworld_runs = []
for filename in filenames:
    import json
    with open(filename, "r") as fp:
        runs = json.load(fp)
    realworld_runs += runs

for run in realworld_runs:
    num_samples = len(run["timestamps"])
    run["timestamps"] = np.concatenate([np.zeros(1), run["timestamps"]])
    run["omega_stator"] = np.concatenate([np.zeros(1), np.ones(num_samples) * run["omega_stator"]])
    run["omega_rotor"] = np.concatenate([np.zeros(1), np.ones(num_samples) * run["omega_rotor"]])
    run["voltages"] = np.concatenate([np.zeros(1), run["Vq"]]) - 1j * np.concatenate([np.zeros(1), run["Vd"]])
    run["currents"] = np.concatenate([np.zeros(1), run["Iq"]]) - 1j * np.concatenate([np.zeros(1), run["Id"]])



# initial guess
motor_params = [1.3, 0.100, 2.0, 0.600, 0.1]

if len(realworld_runs) > 0:
    def test_params(params):
        motor = ACMotor(params[0], params[1], params[2], params[3], params[4])
        residuals = []
        for run in realworld_runs:
            selected_values = run["timestamps"] < 1.0
            
            ys = motor.run(
                run["timestamps"][selected_values],
                run["voltages"].mean(), # TODO: support voltage/omega_s/omega_r series, not just one value
                run["omega_stator"].mean(),
                run["omega_rotor"].mean(),
            )
            
            residuals.append(np.real(ys[1, :] - run["currents"][selected_values]))
            residuals.append(np.imag(ys[1, :] - run["currents"][selected_values]))

        #print(str(params) + " " + str((np.concatenate(residuals)**2).sum()))
        return np.concatenate(residuals)

    result = least_squares(test_params, motor_params)
    if not result.success:
        print("DID NOT CONVERGE!")
    motor_params = result.x



motor = ACMotor(motor_params[0], motor_params[1], motor_params[2], motor_params[3], motor_params[4])

fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
plt.xlabel('Time (s)')
ax1.set_ylabel('Current [A]')
ax2.set_ylabel('Flux [Wb]')

legends = []
for run in realworld_runs:
    selected_values = run["timestamps"] < 2.1
    
    ts = run["timestamps"][selected_values]
    ys = motor.run(
        ts,
        run["voltages"].mean(), # TODO: support voltage/omega_s/omega_r series, not just one value
        run["omega_stator"].mean(),
        run["omega_rotor"].mean(),
    )

    prefix = "$U_s={:.1f}V, \omega_s={:.1f}rad/s$: ".format(run["voltages"].mean(), run["omega_stator"].mean()).replace("+0.0j", "")

    p = ax1.plot(ts, np.real(ys[1, :]), label = prefix + "$I_q$")
    p = ax1.plot(ts, np.real(run["currents"][selected_values]), ".", color=p[0].get_color())
    
    #p = ax1.plot(ts, -np.imag(ys[1, :]), label = prefix + "$I_d$")
    #p = ax1.plot(ts, -np.imag(run["currents"][selected_values]), ".", color=p[0].get_color())

    #ax2.plot(ts, np.absolute(ys[0, :]), prefix + "$\Phi_r$")

fig.legend()
fig.show()


# est. params from using DC dataset:
# [1.23505289, 0.12140287, 2.00802169, 0.55824212, 0.23117988]

# est. params from using blocked rotor dataset:
# [1.56487965, 0.0373343 , 1.95444889, 0.56833555, 0.12249371]

# est. params from using both datasets:
# [1.3948826 , 0.04663297, 1.93612018, 0.58039835, 0.14055246]