From 737808706c21d7030deff04edcdd8dd2f9330321 Mon Sep 17 00:00:00 2001 From: Joel Truher Date: Sun, 15 Oct 2023 10:30:25 -0700 Subject: [PATCH] works with python3, cleaned up the charts a bit --- .gitignore | 2 + control.py | 66 +++++--- dyn_model.py | 433 ++++++++++++++++++++++++-------------------------- misc_utils.py | 31 ++-- my_io.py | 21 +-- my_plot.py | 111 +++++++------ sim_1.py | 128 ++++++++------- 7 files changed, 418 insertions(+), 374 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fffc64d --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +__pycache__ +*.swp diff --git a/control.py b/control.py index f086a7c..d807440 100644 --- a/control.py +++ b/control.py @@ -19,25 +19,26 @@ import numpy as np -import dyn_model as dm +import dyn_model as dm import misc_utils as mu import math PWM_freq = 16000 -PWM_cycle_time = (1./16000) +PWM_cycle_time = 1.0 / 16000 PWM_duty = 0.6 PWM_duty_time = PWM_cycle_time * PWM_duty debug = False + # # # Sp setpoint, Y output # def run_hpwm_l_on_bipol(Sp, Y, t): - elec_angle = mu.norm_angle(Y[dm.ov_theta] * dm.NbPoles/2) + elec_angle = mu.norm_angle(Y[dm.ov_theta] * dm.NbPoles / 2) U = np.zeros(dm.iv_size) @@ -45,7 +46,7 @@ def run_hpwm_l_on_bipol(Sp, Y, t): # switching pattern based on the "encoder" # H PWM L ON pattern - if 0. <= elec_angle <= (math.pi * (1./6.)): # second half of step 1 + if 0.0 <= elec_angle <= (math.pi * (1.0 / 6.0)): # second half of step 1 # U off # V low # W hpwm @@ -59,7 +60,7 @@ def run_hpwm_l_on_bipol(Sp, Y, t): hw = 0 lw = 0 step = "1b" - elif (math.pi * (1.0/6.0)) < elec_angle <= (math.pi * (3.0/6.0)): # step 2 + elif (math.pi * (1.0 / 6.0)) < elec_angle <= (math.pi * (3.0 / 6.0)): # step 2 # U hpwm # V low # W off @@ -73,7 +74,7 @@ def run_hpwm_l_on_bipol(Sp, Y, t): hw = 0 lw = 0 step = "2 " - elif (math.pi * (3.0/6.0)) < elec_angle <= (math.pi * (5.0/6.0)): # step 3 + elif (math.pi * (3.0 / 6.0)) < elec_angle <= (math.pi * (5.0 / 6.0)): # step 3 # U hpwm # V off # W low @@ -87,7 +88,7 @@ def run_hpwm_l_on_bipol(Sp, Y, t): hw = 0 lw = 1 step = "3 " - elif (math.pi * (5.0/6.0)) < elec_angle <= (math.pi * (7.0/6.0)): # step 4 + elif (math.pi * (5.0 / 6.0)) < elec_angle <= (math.pi * (7.0 / 6.0)): # step 4 # U off # V hpwm # W low @@ -101,7 +102,7 @@ def run_hpwm_l_on_bipol(Sp, Y, t): hw = 0 lw = 1 step = "4 " - elif (math.pi * (7.0/6.0)) < elec_angle <= (math.pi * (9.0/6.0)): # step 5 + elif (math.pi * (7.0 / 6.0)) < elec_angle <= (math.pi * (9.0 / 6.0)): # step 5 # U low # V hpwm # W off @@ -115,7 +116,7 @@ def run_hpwm_l_on_bipol(Sp, Y, t): hw = 0 lw = 0 step = "5 " - elif (math.pi * (9.0/6.0)) < elec_angle <= (math.pi * (11.0/6.0)): # step 6 + elif (math.pi * (9.0 / 6.0)) < elec_angle <= (math.pi * (11.0 / 6.0)): # step 6 # U low # V off # W hpwm @@ -129,7 +130,9 @@ def run_hpwm_l_on_bipol(Sp, Y, t): hw = 0 lw = 0 step = "6 " - elif (math.pi * (11.0/6.0)) < elec_angle <= (math.pi * (12.0/6.0)): # first half of step 1 + elif ( + (math.pi * (11.0 / 6.0)) < elec_angle <= (math.pi * (12.0 / 6.0)) + ): # first half of step 1 # U off # V low # W hpwm @@ -144,7 +147,7 @@ def run_hpwm_l_on_bipol(Sp, Y, t): lw = 0 step = "1a" else: - print 'ERROR: The electrical angle is out of range!!!' + print("ERROR: The electrical angle is out of range!!!") # Assigning the scheme phase values to the simulator phases # "Connecting the controller wires to the motor" ^^ @@ -157,16 +160,21 @@ def run_hpwm_l_on_bipol(Sp, Y, t): U[dm.iv_lw] = lv if debug: - print 'time {} step {} eangle {} switches {}'.format(t, step, mu.deg_of_rad(elec_angle), U) + print( + "time {} step {} eangle {} switches {}".format( + t, step, mu.deg_of_rad(elec_angle), U + ) + ) return U + # # # Sp setpoint, Y output # def run_hpwm_l_on(Sp, Y, t): - elec_angle = mu.norm_angle(Y[dm.ov_theta] * dm.NbPoles/2) + elec_angle = mu.norm_angle(Y[dm.ov_theta] * dm.NbPoles / 2) U = np.zeros(dm.iv_size) @@ -174,7 +182,7 @@ def run_hpwm_l_on(Sp, Y, t): # switching pattern based on the "encoder" # H PWM L ON pattern bipolar - if 0. <= elec_angle <= (math.pi * (1./6.)): # second half of step 1 + if 0.0 <= elec_angle <= (math.pi * (1.0 / 6.0)): # second half of step 1 # U off # V low # W hpwm @@ -182,7 +190,7 @@ def run_hpwm_l_on(Sp, Y, t): lu = 0 if math.fmod(t, PWM_cycle_time) <= PWM_duty_time: hw = 1 - lw = 0 + lw = 0 hv = 0 lv = 1 else: @@ -191,24 +199,24 @@ def run_hpwm_l_on(Sp, Y, t): hv = 1 lv = 0 step = "1b" - elif (math.pi * (1.0/6.0)) < elec_angle <= (math.pi * (3.0/6.0)): # step 2 + elif (math.pi * (1.0 / 6.0)) < elec_angle <= (math.pi * (3.0 / 6.0)): # step 2 # U hpwm # V low # W off if math.fmod(t, PWM_cycle_time) <= PWM_duty_time: hu = 1 - lu = 0 + lu = 0 hv = 0 lv = 1 else: hu = 0 - lu = 1 + lu = 1 hv = 1 lv = 0 hw = 0 lw = 0 step = "2 " - elif (math.pi * (3.0/6.0)) < elec_angle <= (math.pi * (5.0/6.0)): # step 3 + elif (math.pi * (3.0 / 6.0)) < elec_angle <= (math.pi * (5.0 / 6.0)): # step 3 # U hpwm # V off # W low @@ -225,7 +233,7 @@ def run_hpwm_l_on(Sp, Y, t): hv = 0 lv = 0 step = "3 " - elif (math.pi * (5.0/6.0)) < elec_angle <= (math.pi * (7.0/6.0)): # step 4 + elif (math.pi * (5.0 / 6.0)) < elec_angle <= (math.pi * (7.0 / 6.0)): # step 4 # U off # V hpwm # W low @@ -242,7 +250,7 @@ def run_hpwm_l_on(Sp, Y, t): hw = 1 lw = 0 step = "4 " - elif (math.pi * (7.0/6.0)) < elec_angle <= (math.pi * (9.0/6.0)): # step 5 + elif (math.pi * (7.0 / 6.0)) < elec_angle <= (math.pi * (9.0 / 6.0)): # step 5 # U low # V hpwm # W off @@ -259,7 +267,7 @@ def run_hpwm_l_on(Sp, Y, t): hw = 0 lw = 0 step = "5 " - elif (math.pi * (9.0/6.0)) < elec_angle <= (math.pi * (11.0/6.0)): # step 6 + elif (math.pi * (9.0 / 6.0)) < elec_angle <= (math.pi * (11.0 / 6.0)): # step 6 # U low # V off # W hpwm @@ -276,7 +284,9 @@ def run_hpwm_l_on(Sp, Y, t): hu = 1 lu = 0 step = "6 " - elif (math.pi * (11.0/6.0)) < elec_angle <= (math.pi * (12.0/6.0)): # first half of step 1 + elif ( + (math.pi * (11.0 / 6.0)) < elec_angle <= (math.pi * (12.0 / 6.0)) + ): # first half of step 1 # U off # V low # W hpwm @@ -294,7 +304,7 @@ def run_hpwm_l_on(Sp, Y, t): lv = 0 step = "1a" else: - print 'ERROR: The electrical angle is out of range!!!' + print("ERROR: The electrical angle is out of range!!!") # Assigning the scheme phase values to the simulator phases # "Connecting the controller wires to the motor" ^^ @@ -307,7 +317,11 @@ def run_hpwm_l_on(Sp, Y, t): U[dm.iv_lw] = lv if debug: - print 'time {} step {} eangle {} switches {}'.format(t, step, mu.deg_of_rad(elec_angle), U) + print( + "time {} step {} eangle {} switches {}".format( + t, step, mu.deg_of_rad(elec_angle), U + ) + ) return U @@ -317,5 +331,5 @@ def run_hpwm_l_on(Sp, Y, t): # Sp setpoint, Y output # def run(Sp, Y, t): - #return run_hpwm_l_on(Sp, Y, t) + # return run_hpwm_l_on(Sp, Y, t) return run_hpwm_l_on_bipol(Sp, Y, t) diff --git a/dyn_model.py b/dyn_model.py index 0a00bfa..fda2de4 100644 --- a/dyn_model.py +++ b/dyn_model.py @@ -1,4 +1,4 @@ -#-*- coding: utf-8 -*- +# -*- coding: utf-8 -*- # # Open-BLDC pysim - Open BrushLess DC Motor Controller python simulator # Copyright (C) 2011 by Antoine Drouin @@ -26,84 +26,84 @@ pset = 2 if pset == 0: - Inertia = 0.0022 # aka. 'J' in kg/(m^2) + Inertia = 0.0022 # aka. 'J' in kg/(m^2) Damping = 0.001 # aka. 'B' in Nm/(rad/s) - Kv = 1700. # aka. motor constant in RPM/V - L = 0.00312 # aka. Coil inductance in H - M = 0.0 # aka. Mutual inductance in H - R = 0.8 # aka. Phase resistence in Ohm - VDC = 100. # aka. Supply voltage - NbPoles = 14. # NbPoles / 2 = Number of pole pairs (you count the permanent magnets on the rotor to get NbPoles) - dvf = .7 # aka. freewheeling diode forward voltage + Kv = 1700.0 # aka. motor constant in RPM/V + L = 0.00312 # aka. Coil inductance in H + M = 0.0 # aka. Mutual inductance in H + R = 0.8 # aka. Phase resistence in Ohm + VDC = 100.0 # aka. Supply voltage + NbPoles = 14.0 # NbPoles / 2 = Number of pole pairs (you count the permanent magnets on the rotor to get NbPoles) + dvf = 0.7 # aka. freewheeling diode forward voltage elif pset == 1: - Inertia = 0.0022 # aka. 'J' in kg/(m^2) + Inertia = 0.0022 # aka. 'J' in kg/(m^2) Damping = 0.001 # aka. 'B' in Nm/(rad/s) - Kv = 70. # aka. motor constant in RPM/V - L = 0.00521 # aka. Coil inductance in H - M = 0.0 # aka. Mutual inductance in H - R = 0.7 # aka. Phase resistence in Ohm - VDC = 100. # aka. Supply voltage - NbPoles = 4. # NbPoles / 2 = Number of pole pairs (you count the permanent magnets on the rotor to get NbPoles) - dvf = .7 # aka. freewheeling diode forward voltage -elif pset == 2: #psim - Inertia = 0.000007 # aka. 'J' in kg/(m^2) + Kv = 70.0 # aka. motor constant in RPM/V + L = 0.00521 # aka. Coil inductance in H + M = 0.0 # aka. Mutual inductance in H + R = 0.7 # aka. Phase resistence in Ohm + VDC = 100.0 # aka. Supply voltage + NbPoles = 4.0 # NbPoles / 2 = Number of pole pairs (you count the permanent magnets on the rotor to get NbPoles) + dvf = 0.7 # aka. freewheeling diode forward voltage +elif pset == 2: # psim + Inertia = 0.000007 # aka. 'J' in kg/(m^2) tau_shaft = 0.006 - Damping = Inertia/tau_shaft # aka. 'B' in Nm/(rad/s) - Kv = 1./32.3*1000 # aka. motor constant in RPM/V - L = 0.00207 # aka. Coil inductance in H - M = -0.00069 # aka. Mutual inductance in H - R = 11.9 # aka. Phase resistence in Ohm - VDC = 100. # aka. Supply voltage - NbPoles = 4. # - dvf = .0 # aka. freewheeling diode forward voltage -elif pset == 3: #modified psim - Inertia = 0.000059 # aka. 'J' in kg/(m^2) + Damping = Inertia / tau_shaft # aka. 'B' in Nm/(rad/s) + Kv = 1.0 / 32.3 * 1000 # aka. motor constant in RPM/V + L = 0.00207 # aka. Coil inductance in H + M = -0.00069 # aka. Mutual inductance in H + R = 11.9 # aka. Phase resistence in Ohm + VDC = 100.0 # aka. Supply voltage + NbPoles = 4.0 # + dvf = 0.0 # aka. freewheeling diode forward voltage +elif pset == 3: # modified psim + Inertia = 0.000059 # aka. 'J' in kg/(m^2) tau_shaft = 0.006 - Damping = Inertia/tau_shaft # aka. 'B' in Nm/(rad/s) - Kv = 1./32.3*1000 # aka. motor constant in RPM/V - L = 0.00207 # aka. Coil inductance in H - M = -0.00069 # aka. Mutual inductance in H - R = 11.9 # aka. Phase resistence in Ohm - VDC = 300. # aka. Supply voltage - NbPoles = 4. # - dvf = .0 # aka. freewheeling diode forward voltage + Damping = Inertia / tau_shaft # aka. 'B' in Nm/(rad/s) + Kv = 1.0 / 32.3 * 1000 # aka. motor constant in RPM/V + L = 0.00207 # aka. Coil inductance in H + M = -0.00069 # aka. Mutual inductance in H + R = 11.9 # aka. Phase resistence in Ohm + VDC = 300.0 # aka. Supply voltage + NbPoles = 4.0 # + dvf = 0.0 # aka. freewheeling diode forward voltage else: - print "Unknown pset {}".format(pset) + print("Unknown pset {}".format(pset)) # Components of the state vector -sv_theta = 0 # angle of the rotor -sv_omega = 1 # angular speed of the rotor -sv_iu = 2 # phase u current -sv_iv = 3 # phase v current -sv_iw = 4 # phase w current +sv_theta = 0 # angle of the rotor +sv_omega = 1 # angular speed of the rotor +sv_iu = 2 # phase u current +sv_iv = 3 # phase v current +sv_iw = 4 # phase w current sv_size = 5 # Components of the command vector -iv_lu = 0 -iv_hu = 1 -iv_lv = 2 -iv_hv = 3 -iv_lw = 4 -iv_hw = 5 +iv_lu = 0 +iv_hu = 1 +iv_lv = 2 +iv_hv = 3 +iv_lw = 4 +iv_hw = 5 iv_size = 6 # Components of the perturbation vector -pv_torque = 0 +pv_torque = 0 pv_friction = 1 pv_size = 2 # Components of the output vector -ov_iu = 0 -ov_iv = 1 -ov_iw = 2 -ov_vu = 3 -ov_vv = 4 -ov_vw = 5 -ov_theta = 6 -ov_omega = 7 -ov_size = 8 +ov_iu = 0 +ov_iv = 1 +ov_iw = 2 +ov_vu = 3 +ov_vv = 4 +ov_vw = 5 +ov_theta = 6 +ov_omega = 7 +ov_size = 8 # Phases and star vector designators ph_U = 0 @@ -122,174 +122,173 @@ dv_ph_star = 6 dv_size = 7 + # # Calculate backemf at a given omega offset from the current rotor position # # Used to calculate the phase backemf aka. 'e' # -def backemf(X,thetae_offset): - phase_thetae = mu.norm_angle((X[sv_theta] * (NbPoles / 2.)) + thetae_offset) +def backemf(X, thetae_offset): + phase_thetae = mu.norm_angle((X[sv_theta] * (NbPoles / 2.0)) + thetae_offset) - bemf_constant = mu.vpradps_of_rpmpv(Kv) # aka. ke in V/rad/s + bemf_constant = mu.vpradps_of_rpmpv(Kv) # aka. ke in V/rad/s max_bemf = bemf_constant * X[sv_omega] - bemf = 0. - if 0. <= phase_thetae <= (math.pi * (1./6.)): - bemf = (max_bemf / (math.pi * (1./6.))) * phase_thetae - elif (math.pi/6.) < phase_thetae <= (math.pi * (5./6.)): + bemf = 0.0 + if 0.0 <= phase_thetae <= (math.pi * (1.0 / 6.0)): + bemf = (max_bemf / (math.pi * (1.0 / 6.0))) * phase_thetae + elif (math.pi / 6.0) < phase_thetae <= (math.pi * (5.0 / 6.0)): bemf = max_bemf - elif (math.pi * (5./6.)) < phase_thetae <= (math.pi * (7./6.)): - bemf = -((max_bemf/(math.pi/6.))* (phase_thetae - math.pi)) - elif (math.pi * (7./6.)) < phase_thetae <= (math.pi * (11./6.)): + elif (math.pi * (5.0 / 6.0)) < phase_thetae <= (math.pi * (7.0 / 6.0)): + bemf = -((max_bemf / (math.pi / 6.0)) * (phase_thetae - math.pi)) + elif (math.pi * (7.0 / 6.0)) < phase_thetae <= (math.pi * (11.0 / 6.0)): bemf = -max_bemf - elif (math.pi * (11./6.)) < phase_thetae <= (2.0 * math.pi): - bemf = (max_bemf/(math.pi/6.)) * (phase_thetae - (2. * math.pi)) + elif (math.pi * (11.0 / 6.0)) < phase_thetae <= (2.0 * math.pi): + bemf = (max_bemf / (math.pi / 6.0)) * (phase_thetae - (2.0 * math.pi)) else: - print "ERROR: angle out of bounds can not calculate bemf {}".format(phase_thetae) + print( + "ERROR: angle out of bounds can not calculate bemf {}".format(phase_thetae) + ) return bemf + # # Calculate phase voltages # Returns a vector of phase voltages in reference to the star point def voltages(X, U): - - eu = backemf(X, 0.) - ev = backemf(X, math.pi * (2./3.)) - ew = backemf(X, math.pi * (4./3.)) + eu = backemf(X, 0.0) + ev = backemf(X, math.pi * (2.0 / 3.0)) + ew = backemf(X, math.pi * (4.0 / 3.0)) # Check which phases are excited - pux = (U[iv_hu] == 1) or \ - (U[iv_lu] == 1) + pux = (U[iv_hu] == 1) or (U[iv_lu] == 1) - pvx = (U[iv_hv] == 1) or \ - (U[iv_lv] == 1) + pvx = (U[iv_hv] == 1) or (U[iv_lv] == 1) - pwx = (U[iv_hw] == 1) or \ - (U[iv_lw] == 1) + pwx = (U[iv_hw] == 1) or (U[iv_lw] == 1) - vu = 0. - vv = 0. - vw = 0. - vm = 0. + vu = 0.0 + vv = 0.0 + vw = 0.0 + vm = 0.0 if pux and pvx and pwx: - if (U[iv_hu] == 1): - vu = VDC/2. + if U[iv_hu] == 1: + vu = VDC / 2.0 else: - vu = -VDC/2. + vu = -VDC / 2.0 - if (U[iv_hv] == 1): - vv = VDC/2. + if U[iv_hv] == 1: + vv = VDC / 2.0 else: - vv = -VDC/2. + vv = -VDC / 2.0 - if (U[iv_hw] == 1): - vw = VDC/2. + if U[iv_hw] == 1: + vw = VDC / 2.0 else: - vw = -VDC/2. + vw = -VDC / 2.0 - vm = (vu + vv + vw - eu - ev - ew) / 3. + vm = (vu + vv + vw - eu - ev - ew) / 3.0 elif pux and pvx: - # calculate excited phase voltages - if (U[iv_hu] == 1): - vu = VDC/2. + if U[iv_hu] == 1: + vu = VDC / 2.0 else: - vu = -VDC/2. + vu = -VDC / 2.0 - if (U[iv_hv] == 1): - vv = VDC/2. + if U[iv_hv] == 1: + vv = VDC / 2.0 else: - vv = -VDC/2. + vv = -VDC / 2.0 # calculate star voltage - vm = (vu + vv - eu - ev) / 2. + vm = (vu + vv - eu - ev) / 2.0 # calculate remaining phase voltage vw = ew + vm # clip the voltage to freewheeling diodes - #if (vw > ((VDC/2) + dvf)): + # if (vw > ((VDC/2) + dvf)): # vw = (VDC/2) + dvf; # vm = (vu + vv + vw - eu - ev - ew) / 3. - #elif (vw < (-(VDC/2) - dvf)): + # elif (vw < (-(VDC/2) - dvf)): # vw = -(VDC/2) - dvf; # vm = (vu + vv + vw - eu - ev - ew) / 3. elif pux and pwx: - if (U[iv_hu] == 1): - vu = VDC/2. + if U[iv_hu] == 1: + vu = VDC / 2.0 else: - vu = -VDC/2. + vu = -VDC / 2.0 - if (U[iv_hw] == 1): - vw = VDC/2. + if U[iv_hw] == 1: + vw = VDC / 2.0 else: - vw = -VDC/2. + vw = -VDC / 2.0 - vm = (vu + vw - eu - ew) / 2. + vm = (vu + vw - eu - ew) / 2.0 vv = ev + vm # clip the voltage to freewheeling diodes - #if (vv > ((VDC/2) + dvf)): + # if (vv > ((VDC/2) + dvf)): # vv = (VDC/2) + dvf; # vm = (vu + vv + vw - eu - ev - ew) / 3. - #elif (vv < (-(VDC/2) - dvf)): + # elif (vv < (-(VDC/2) - dvf)): # vv = -(VDC/2) - dvf; # vm = (vu + vv + vw - eu - ev - ew) / 3. elif pvx and pwx: - if (U[iv_hv] == 1): - vv = VDC/2. + if U[iv_hv] == 1: + vv = VDC / 2.0 else: - vv = -VDC/2. + vv = -VDC / 2.0 - if (U[iv_hw] == 1): - vw = VDC/2. + if U[iv_hw] == 1: + vw = VDC / 2.0 else: - vw = -VDC/2. + vw = -VDC / 2.0 - vm = (vv + vw - ev - ew) / 2. + vm = (vv + vw - ev - ew) / 2.0 vu = eu + vm # clip the voltage to freewheeling diodes - #if (vu > ((VDC/2) + dvf)): + # if (vu > ((VDC/2) + dvf)): # vu = (VDC/2) + dvf; # vm = (vu + vv + vw - eu - ev - ew) / 3. - #elif (vu < (-(VDC/2) - dvf)): + # elif (vu < (-(VDC/2) - dvf)): # vu = -(VDC/2) - dvf; # vm = (vu + vv + vw - eu - ev - ew) / 3. elif pux: - if (U[iv_hu] == 1): - vu = VDC/2 + if U[iv_hu] == 1: + vu = VDC / 2 else: - vu = -VDC/2. + vu = -VDC / 2.0 - vm = (vu - eu) + vm = vu - eu vv = ev + vm vw = ew + vm - # if we want to handle diodes properly how to do that here? + # if we want to handle diodes properly how to do that here? elif pvx: - if (U[iv_hv] == 1): - vv = VDC/2 + if U[iv_hv] == 1: + vv = VDC / 2 else: - vv = -VDC/2. + vv = -VDC / 2.0 - vm = (vv - ev) + vm = vv - ev vu = eu + vm vw = ew + vm elif pwx: - if (U[iv_hw] == 1): - vw = VDC/2 + if U[iv_hw] == 1: + vw = VDC / 2 else: - vw = -VDC/2. + vw = -VDC / 2.0 - vm = (vw - ew) + vm = vw - ew vu = eu + vm vv = ev + vm else: @@ -297,62 +296,57 @@ def voltages(X, U): vv = ev vw = ew - -# # Initialize the imposed terminal voltages -# vui = 0. -# vvi = 0. -# vwi = 0. -# -# # Phase input voltages based on the inverter switches states -# if (U[iv_hu] == 1) or (U[iv_dhu] == 1): -# vui = VDC/2. -# if (U[iv_lu] == 1) or (U[iv_dlu] == 1): -# vui = -VDC/2. -# if (U[iv_hv] == 1) or (U[iv_dhv] == 1): -# vvi = VDC/2. -# if (U[iv_lv] == 1) or (U[iv_dlv] == 1): -# vvi = -VDC/2. -# if (U[iv_hw] == 1) or (U[iv_dhw] == 1): -# vwi = VDC/2. -# if (U[iv_lw] == 1) or (U[iv_dlw] == 1): -# vwi = -VDC/2. -# -# #i_thr = 0.001 # current threshold saying that the phase is not conducting -# i_thr = 0. # current threshold saying that the phase is not conducting -# #if -i_thr < X[sv_iu] < i_thr: # phase V & W are conducting current -# if not pux: # phase V & W are conducting current -# vm = ((vvi + vwi) / 2.) - ((ev + ew) / 2.) -# vu = eu -# vv = vvi - vm -# vw = vwi - vm -# elif not pvx: # phase U & W are conducting current -# vm = ((vui + vwi) / 2.) - ((eu + ew) / 2.) -# vu = vui - vm -# vv = ev -# vw = vwi - vm -# elif not pwx: # phase U & V are conducting current -# vm = ((vui + vvi) / 2.) - ((eu + ev) / 2.) -# vu = vui - vm -# vv = vvi - vm -# vw = ew -# else: # all phases are corducting current -# print "all phases are conducting!" -# vm = ((vui + vvi + vwi) / 3.) - ((eu + ev + ew) / 3.) -# vu = vui - vm -# vv = vvi - vm -# vw = vwi - vm - - -# print "{} : {} {} {}".format(X[sv_omega], vu, vv, vw ) - - V = [ vu, - vv, - vw, - vm - ] + # # Initialize the imposed terminal voltages + # vui = 0. + # vvi = 0. + # vwi = 0. + # + # # Phase input voltages based on the inverter switches states + # if (U[iv_hu] == 1) or (U[iv_dhu] == 1): + # vui = VDC/2. + # if (U[iv_lu] == 1) or (U[iv_dlu] == 1): + # vui = -VDC/2. + # if (U[iv_hv] == 1) or (U[iv_dhv] == 1): + # vvi = VDC/2. + # if (U[iv_lv] == 1) or (U[iv_dlv] == 1): + # vvi = -VDC/2. + # if (U[iv_hw] == 1) or (U[iv_dhw] == 1): + # vwi = VDC/2. + # if (U[iv_lw] == 1) or (U[iv_dlw] == 1): + # vwi = -VDC/2. + # + # #i_thr = 0.001 # current threshold saying that the phase is not conducting + # i_thr = 0. # current threshold saying that the phase is not conducting + # #if -i_thr < X[sv_iu] < i_thr: # phase V & W are conducting current + # if not pux: # phase V & W are conducting current + # vm = ((vvi + vwi) / 2.) - ((ev + ew) / 2.) + # vu = eu + # vv = vvi - vm + # vw = vwi - vm + # elif not pvx: # phase U & W are conducting current + # vm = ((vui + vwi) / 2.) - ((eu + ew) / 2.) + # vu = vui - vm + # vv = ev + # vw = vwi - vm + # elif not pwx: # phase U & V are conducting current + # vm = ((vui + vvi) / 2.) - ((eu + ev) / 2.) + # vu = vui - vm + # vv = vvi - vm + # vw = ew + # else: # all phases are corducting current + # print "all phases are conducting!" + # vm = ((vui + vvi + vwi) / 3.) - ((eu + ev + ew) / 3.) + # vu = vui - vm + # vv = vvi - vm + # vw = vwi - vm + + # print "{} : {} {} {}".format(X[sv_omega], vu, vv, vw ) + + V = [vu, vv, vw, vm] return V + # # Dynamic model # @@ -363,55 +357,42 @@ def dyn(X, t, U, W): return Xd + # Dynamic model with debug vector def dyn_debug(X, t, U, W): - - eu = backemf(X, 0.) - ev = backemf(X, math.pi * (2./3.)) - ew = backemf(X, math.pi * (4./3.)) + eu = backemf(X, 0.0) + ev = backemf(X, math.pi * (2.0 / 3.0)) + ew = backemf(X, math.pi * (4.0 / 3.0)) # Electromagnetic torque - etorque = (eu * X[sv_iu] + ev * X[sv_iv] + ew * X[sv_iw])/X[sv_omega] + etorque = (eu * X[sv_iu] + ev * X[sv_iv] + ew * X[sv_iw]) / X[sv_omega] # Mechanical torque - mtorque = ((etorque * (NbPoles / 2)) - (Damping * X[sv_omega]) - W[pv_torque]) + mtorque = (etorque * (NbPoles / 2)) - (Damping * X[sv_omega]) - W[pv_torque] - if ((mtorque > 0) and (mtorque <= W[pv_friction])): + if (mtorque > 0) and (mtorque <= W[pv_friction]): mtorque = 0 - elif (mtorque >= W[pv_friction]): + elif mtorque >= W[pv_friction]: mtorque = mtorque - W[pv_friction] - elif ((mtorque < 0) and (mtorque >= (-W[pv_friction]))): - mtorque = 0 - elif (mtorque <= (-W[pv_friction])): - mtorque = mtorque + W[pv_friction] + elif (mtorque < 0) and (mtorque >= (-W[pv_friction])): + mtorque = 0 + elif mtorque <= (-W[pv_friction]): + mtorque = mtorque + W[pv_friction] # Acceleration of the rotor omega_dot = mtorque / Inertia V = voltages(X, U) - pdt = VDC/2 + dvf + pdt = VDC / 2 + dvf iu_dot = (V[ph_U] - (R * X[sv_iu]) - eu - V[ph_star]) / (L - M) iv_dot = (V[ph_V] - (R * X[sv_iv]) - ev - V[ph_star]) / (L - M) iw_dot = (V[ph_W] - (R * X[sv_iw]) - ew - V[ph_star]) / (L - M) - Xd = [ X[sv_omega], - omega_dot, - iu_dot, - iv_dot, - iw_dot - ] - - Xdebug = [ - eu, - ev, - ew, - V[ph_U], - V[ph_V], - V[ph_W], - V[ph_star] - ] + Xd = [X[sv_omega], omega_dot, iu_dot, iv_dot, iw_dot] + + Xdebug = [eu, ev, ew, V[ph_U], V[ph_V], V[ph_W], V[ph_star]] return Xd, Xdebug @@ -420,11 +401,17 @@ def dyn_debug(X, t, U, W): # # def output(X, U): - V = voltages(X, U) - Y = [X[sv_iu], X[sv_iv], X[sv_iw], - V[ph_U], V[ph_V], V[ph_W], - X[sv_theta], X[sv_omega]] + Y = [ + X[sv_iu], + X[sv_iv], + X[sv_iw], + V[ph_U], + V[ph_V], + V[ph_W], + X[sv_theta], + X[sv_omega], + ] return Y diff --git a/misc_utils.py b/misc_utils.py index e3025ff..f857f32 100644 --- a/misc_utils.py +++ b/misc_utils.py @@ -19,29 +19,42 @@ import math + # -def rad_of_deg(d): return d/180.*math.pi +def rad_of_deg(d): + return d / 180.0 * math.pi + # -def deg_of_rad(r): return r*180./math.pi +def deg_of_rad(r): + return r * 180.0 / math.pi + # -def rpm_of_radps(rps): return rps/(2*math.pi)*60 +def rpm_of_radps(rps): + return rps / (2 * math.pi) * 60 + # -def degps_of_radps(rps): return rps/(2*math.pi)*60*360 +def degps_of_radps(rps): + return rps / (2 * math.pi) * 60 * 360 + # -def radps_of_rpm(rpm): return rpm*(2*math.pi)/60 +def radps_of_rpm(rpm): + return rpm * (2 * math.pi) / 60 + # -def vpradps_of_rpmpv(vprpm): return 30/(vprpm*math.pi) +def vpradps_of_rpmpv(vprpm): + return 30 / (vprpm * math.pi) + # def norm_angle(alpha): - alpha_n = math.fmod(alpha, 2*math.pi) + alpha_n = math.fmod(alpha, 2 * math.pi) - if alpha_n < 0.: - alpha_n = (2*math.pi) + alpha_n + if alpha_n < 0.0: + alpha_n = (2 * math.pi) + alpha_n return alpha_n diff --git a/my_io.py b/my_io.py index 12917bd..915d86a 100644 --- a/my_io.py +++ b/my_io.py @@ -17,19 +17,20 @@ # along with this program. If not, see . # -import numpy as np -import dyn_model as dm +import numpy as np +import dyn_model as dm import misc_utils as mu + def read_csv(filename): my_records = np.recfromcsv(filename) Y = np.zeros((my_records.time.size, dm.ov_size)) - Y[:,dm.ov_iu] = my_records.ia - Y[:,dm.ov_iv] = my_records.ib - Y[:,dm.ov_iw] = my_records.ic - Y[:,dm.ov_vu] = my_records.vag - Y[:,dm.ov_vv] = my_records.vbg - Y[:,dm.ov_vw] = my_records.vcg - Y[:,dm.ov_omega] = mu.radps_of_rpm(my_records.nm) -# import code; code.interact(local=locals()) + Y[:, dm.ov_iu] = my_records.ia + Y[:, dm.ov_iv] = my_records.ib + Y[:, dm.ov_iw] = my_records.ic + Y[:, dm.ov_vu] = my_records.vag + Y[:, dm.ov_vv] = my_records.vbg + Y[:, dm.ov_vw] = my_records.vcg + Y[:, dm.ov_omega] = mu.radps_of_rpm(my_records.nm) + # import code; code.interact(local=locals()) return my_records.time, Y diff --git a/my_plot.py b/my_plot.py index 1e98e75..9b13102 100644 --- a/my_plot.py +++ b/my_plot.py @@ -19,81 +19,96 @@ import matplotlib.pyplot as plt -import dyn_model as dm +import dyn_model as dm import misc_utils as mu ang_unit_rad_s = 0 ang_unit_deg_s = 1 ang_unit_rpm = 2 + def plot_output(time, Y, ls): ang_unit = ang_unit_rpm + fig, ax = plt.subplots(4, 1, figsize=(8,8)) + fig.suptitle('Output') + # Phase current - ax = plt.subplot(4, 1, 1) - ax.yaxis.set_label_text('A', {'color' : 'k', 'fontsize' : 15 }) - plt.plot(time,Y[:,dm.ov_iu], ls, linewidth=1.5) - plt.plot(time,Y[:,dm.ov_iv], ls, linewidth=1.5) - plt.plot(time,Y[:,dm.ov_iw], ls, linewidth=1.5) - plt.legend(['$i_u$', '$i_v$', '$i_w$'], loc='upper right') - plt.title('Phase current') + #ax = plt.subplot(4, 1, 1) + ax[0].yaxis.set_label_text("A", {"color": "k", "fontsize": 15}) + ax[0].plot(time, Y[:, dm.ov_iu], ls, linewidth=1.5) + ax[0].plot(time, Y[:, dm.ov_iv], ls, linewidth=1.5) + ax[0].plot(time, Y[:, dm.ov_iw], ls, linewidth=1.5) + ax[0].legend(["$i_u$", "$i_v$", "$i_w$"], loc="upper right") + ax[0].set_title("Phase current") # Phase terminal voltage - ax = plt.subplot(4, 1, 2) - ax.yaxis.set_label_text('V', {'color' : 'k', 'fontsize' : 15 }) - plt.plot(time,Y[:,dm.ov_vu], ls, linewidth=1.5) - plt.plot(time,Y[:,dm.ov_vv], ls, linewidth=1.5) - plt.plot(time,Y[:,dm.ov_vw], ls, linewidth=1.5) - plt.legend(['$v_u$', '$v_v$', '$v_w$'], loc='upper right') - plt.title('Phase terminal voltage') + #ax = plt.subplot(4, 1, 2) + ax[1].yaxis.set_label_text("V", {"color": "k", "fontsize": 15}) + ax[1].plot(time, Y[:, dm.ov_vu], ls, linewidth=1.5) + ax[1].plot(time, Y[:, dm.ov_vv], ls, linewidth=1.5) + ax[1].plot(time, Y[:, dm.ov_vw], ls, linewidth=1.5) + ax[1].legend(["$v_u$", "$v_v$", "$v_w$"], loc="upper right") + ax[1].set_title("Phase terminal voltage") # Rotor mechanical position - ax = plt.subplot(4, 1, 3) - ax.yaxis.set_label_text('Deg', {'color' : 'k', 'fontsize' : 15 }) - plt.plot(time,mu.deg_of_rad(Y[:,dm.ov_theta]), ls, linewidth=1.5) -# plt.plot(time, Y[:,dm.ov_theta], ls, linewidth=1.5) - plt.title('Rotor angular position') + #ax = plt.subplot(4, 1, 3) + ax[2].yaxis.set_label_text("Deg", {"color": "k", "fontsize": 15}) + ax[2].plot(time, mu.deg_of_rad(Y[:, dm.ov_theta]), ls, linewidth=1.5) + # plt.plot(time, Y[:,dm.ov_theta], ls, linewidth=1.5) + ax[2].set_title("Rotor angular position") # Rotor mechanical angular speed - ax = plt.subplot(4, 1, 4) + #ax = plt.subplot(4, 1, 4) + + if ang_unit == ang_unit_rad_s: + ax[3].yaxis.set_label_text("Rad/s", {"color": "k", "fontsize": 15}) + ax[3].plot(time, Y[:, dm.ov_omega], ls, linewidth=1.5) + elif ang_unit == ang_unit_deg_s: + ax[3].yaxis.set_label_text("Deg/s", {"color": "k", "fontsize": 15}) + ax[3].plot(time, mu.degps_of_radps(Y[:, dm.ov_omega]), ls, linewidth=1.5) + elif ang_unit == ang_unit_rpm: + ax[3].yaxis.set_label_text("RPM", {"color": "k", "fontsize": 15}) + ax[3].plot(time, mu.rpm_of_radps(Y[:, dm.ov_omega]), ls, linewidth=1.5) + + ax[3].set_title("Rotor Rotational Velocity") - if (ang_unit == ang_unit_rad_s): - ax.yaxis.set_label_text('Rad/s', {'color' : 'k', 'fontsize' : 15 }) - plt.plot(time,Y[:,dm.ov_omega], ls, linewidth=1.5) - elif (ang_unit == ang_unit_deg_s): - ax.yaxis.set_label_text('Deg/s', {'color' : 'k', 'fontsize' : 15 }) - plt.plot(time,mu.degps_of_radps(Y[:,dm.ov_omega]), ls, linewidth=1.5) - elif (ang_unit == ang_unit_rpm): - ax.yaxis.set_label_text('RPM', {'color' : 'k', 'fontsize' : 15 }) - plt.plot(time,mu.rpm_of_radps(Y[:,dm.ov_omega]), ls, linewidth=1.5) + fig.tight_layout() - plt.title('Rotor Rotational Velocity') def plot_debug(time, Xdebug): - plt.subplot(4, 1, 1) + fig, ax = plt.subplots(3, 1, figsize=(8,8)) + fig.suptitle('Debug') + + #plt.subplot(4, 1, 1) + + ax[0].plot(time, Xdebug[:, dm.dv_eu], linewidth=1.5) + ax[0].plot(time, Xdebug[:, dm.dv_ev], linewidth=1.5) + ax[0].plot(time, Xdebug[:, dm.dv_ew], linewidth=1.5) + ax[0].legend(["$U_{BEMF}$", "$V_{BEMF}$", "$W_{BEMF}$"], loc="upper right") + + #plt.subplot(4, 1, 2) - plt.plot(time,Xdebug[:,dm.dv_eu], linewidth=1.5) - plt.plot(time,Xdebug[:,dm.dv_ev], linewidth=1.5) - plt.plot(time,Xdebug[:,dm.dv_ew], linewidth=1.5) - plt.legend(['$U_{BEMF}$', '$V_{BEMF}$', '$W_{BEMF}$'], loc='upper right') + ax[1].plot(time, Xdebug[:, dm.dv_ph_U], linewidth=1.5) + ax[1].plot(time, Xdebug[:, dm.dv_ph_V], linewidth=1.5) + ax[1].plot(time, Xdebug[:, dm.dv_ph_W], linewidth=1.5) + ax[1].legend(["$U$", "$V$", "$W$"], loc="upper right") - plt.subplot(4, 1, 2) + #plt.subplot(4, 1, 3) - plt.plot(time,Xdebug[:,dm.dv_ph_U], linewidth=1.5) - plt.plot(time,Xdebug[:,dm.dv_ph_V], linewidth=1.5) - plt.plot(time,Xdebug[:,dm.dv_ph_W], linewidth=1.5) - plt.legend(['$U$', '$V$', '$W$'], loc='upper right') + ax[2].plot(time, Xdebug[:, dm.dv_ph_star], linewidth=1.5) + ax[2].legend(["$star$"], loc="upper right") - plt.subplot(4, 1, 3) + fig.tight_layout() - plt.plot(time,Xdebug[:,dm.dv_ph_star], linewidth=1.5) - plt.legend(['$star$'], loc='upper right') def plot_diodes(time, D): + titles_diodes = ["$dhu$", "$dlu$", "$dhv$", "$dlv$", "$dhw$", "$dlw$"] - titles_diodes = ['$dhu$', '$dlu$', '$dhv$', '$dlv$', '$dhw$', '$dlw$'] + fig, ax = plt.subplots(6, 2) + fig.suptitle('Diodes') for i in range(0, dm.adc_size): - plt.subplot(6, 2, 2*i+1) - plt.plot(time,D[:,i], 'r', linewidth=1.5) - plt.title(titles_diodes[i]) + #plt.subplot(6, 2, 2 * i + 1) + ax[i][0].plot(time, D[:, i], "r", linewidth=1.5) + ax[i][0].title(titles_diodes[i]) diff --git a/sim_1.py b/sim_1.py index 12ef53a..cdc5e44 100755 --- a/sim_1.py +++ b/sim_1.py @@ -1,4 +1,3 @@ - # # Open-BLDC pysim - Open BrushLess DC Motor Controller python simulator # Copyright (C) 2011 by Antoine Drouin @@ -19,87 +18,101 @@ # import matplotlib -matplotlib.use('MacOSX') -#matplotlib.use('GTKCairo') import numpy as np -import pylab as pl import matplotlib.pyplot as plt from scipy import integrate from scipy.signal import decimate import misc_utils as mu -import dyn_model as dm -import control as ctl -import my_io as mio -import my_plot as mp - +import dyn_model as dm +import control as ctl +import my_io as mio +import my_plot as mp def display_state_and_command(time, X, U): + titles_state = ["$\\theta$", "$\\omega$", "$i_u$", "$i_v$", "$i_w$"] + titles_cmd = ["$u_l$", "$u_h$", "$v_l$", "$v_h$", "$w_l$", "$w_h$"] + + #plt.figure(figsize=(10.24, 5.12)) + fig, ax = plt.subplots(6, 2, figsize=(10,10)) + fig.suptitle('State and Command') - titles_state = ['$\\theta$', '$\omega$', '$i_u$', '$i_v$', '$i_w$'] - titles_cmd = ['$u_l$', '$u_h$', '$v_l$', '$v_h$', '$w_l$', '$w_h$'] for i in range(0, 2): - plt.subplot(6, 2, 2*i+1) - plt.plot(time,mu.deg_of_rad(X[:,i]), 'r', linewidth=3.0) - plt.title(titles_state[i]) + #plt.subplot(6, 2, 2 * i + 1) + ax[i][0].plot(time, mu.deg_of_rad(X[:, i]), "r", linewidth=3.0) + ax[i][0].set_title(titles_state[i]) for i in range(2, dm.sv_size): - plt.subplot(6, 2, 2*i+1) - plt.plot(time, X[:,i], 'r', linewidth=3.0) - plt.title(titles_state[i]) + #plt.subplot(6, 2, 2 * i + 1) + ax[i][0].plot(time, X[:, i], "r", linewidth=3.0) + ax[i][0].set_title(titles_state[i]) for i in range(0, 6): - plt.subplot(6, 2, 2*i+2) - plt.plot(time, U[:,i], 'r', linewidth=3.0) - plt.title(titles_cmd[i]) + #plt.subplot(6, 2, 2 * i + 2) + ax[i][1].plot(time, U[:, i], "r", linewidth=3.0) + ax[i][1].set_title(titles_cmd[i]) + + fig.tight_layout() + + def print_simulation_progress(count, steps): - sim_perc_last = ((count-1)*100) / steps - sim_perc = (count*100) / steps - if (sim_perc_last != sim_perc): - print "{}%".format(sim_perc) + sim_perc_last = ((count - 1) * 100) / steps + sim_perc = (count * 100) / steps + if sim_perc_last != sim_perc: + print("{}%".format(sim_perc)) + def drop_it(a, factor): - new = [] - for n, x in enumerate(a): - if ((n % factor) == 0): - new.append(x) - return np.array(new) + new = [] + for n, x in enumerate(a): + if (n % factor) == 0: + new.append(x) + return np.array(new) + def compress(a, factor): - return drop_it(a, factor) - #return decimate(a, 8, n=8, axis=0) + return drop_it(a, factor) + # return decimate(a, 8, n=8, axis=0) + def main(): -# t_psim, Y_psim = mio.read_csv('bldc_startup_psim_1us_resolution.csv') -# mp.plot_output(t_psim, Y_psim, '.') + # t_psim, Y_psim = mio.read_csv('bldc_startup_psim_1us_resolution.csv') + # mp.plot_output(t_psim, Y_psim, '.') - freq_sim = 1e6 # simulation frequency + freq_sim = 1e6 # simulation frequency compress_factor = 3 - time = pl.arange(0.0, 0.01, 1./freq_sim) # create time slice vector - X = np.zeros((time.size, dm.sv_size)) # allocate state vector + time = np.arange(0.0, 0.01, 1.0 / freq_sim) # create time slice vector + X = np.zeros((time.size, dm.sv_size)) # allocate state vector Xdebug = np.zeros((time.size, dm.dv_size)) # allocate debug data vector - Y = np.zeros((time.size, dm.ov_size)) # allocate output vector - U = np.zeros((time.size, dm.iv_size)) # allocate input vector - X0 = [0, mu.rad_of_deg(0.1), 0, 0, 0] # - X[0,:] = X0 + Y = np.zeros((time.size, dm.ov_size)) # allocate output vector + U = np.zeros((time.size, dm.iv_size)) # allocate input vector + X0 = [0, mu.rad_of_deg(0.1), 0, 0, 0] # + X[0, :] = X0 W = [0, 1] - for i in range(1,time.size): - - if i==1: + for i in range(1, time.size): + if i == 1: Uim2 = np.zeros(dm.iv_size) else: - Uim2 = U[i-2,:] - - Y[i-1,:] = dm.output(X[i-1,:], Uim2) # get the output for the last step - U[i-1,:] = ctl.run(0, Y[i-1,:], time[i-1]) # run the controller for the last step - tmp = integrate.odeint(dm.dyn, X[i-1,:], [time[i-1], time[i]], args=(U[i-1,:], W)) # integrate - X[i,:] = tmp[1,:] # copy integration output to the current step - X[i, dm.sv_theta] = mu.norm_angle( X[i, dm.sv_theta]) # normalize the angle in the state - tmp, Xdebug[i,:] = dm.dyn_debug(X[i-1,:], time[i-1], U[i-1,:], W) # get debug data + Uim2 = U[i - 2, :] + + Y[i - 1, :] = dm.output(X[i - 1, :], Uim2) # get the output for the last step + U[i - 1, :] = ctl.run( + 0, Y[i - 1, :], time[i - 1] + ) # run the controller for the last step + tmp = integrate.odeint( + dm.dyn, X[i - 1, :], [time[i - 1], time[i]], args=(U[i - 1, :], W) + ) # integrate + X[i, :] = tmp[1, :] # copy integration output to the current step + X[i, dm.sv_theta] = mu.norm_angle( + X[i, dm.sv_theta] + ) # normalize the angle in the state + tmp, Xdebug[i, :] = dm.dyn_debug( + X[i - 1, :], time[i - 1], U[i - 1, :], W + ) # get debug data print_simulation_progress(i, time.size) - Y[-1,:] = Y[-2,:] - U[-1,:] = U[-2,:] + Y[-1, :] = Y[-2, :] + U[-1, :] = U[-2, :] if compress_factor > 1: time = compress(time, compress_factor) @@ -108,15 +121,14 @@ def main(): U = compress(U, compress_factor) Xdebug = compress(Xdebug, compress_factor) - mp.plot_output(time, Y, '-') -# pl.show() - plt.figure(figsize=(10.24, 5.12)) + mp.plot_output(time, Y, "-") display_state_and_command(time, X, U) - plt.figure(figsize=(10.24, 5.12)) + #plt.figure(figsize=(10.24, 5.12)) mp.plot_debug(time, Xdebug) - pl.show() + plt.show() + if __name__ == "__main__": main()