DriCy_Selo.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Created on Thu Mar  4 08:53:28 2021

@author:    Michael Schlüter \n
            Technische Universität Berlin \n
            https://www.pe.tu-berlin.de/ \n
            m.schlueter@tu-berlin.de
"""


from unittest.loader import VALID_MODULE_NAME
import toml #used for conifg data (vehicles, semiconductor, ...)
import pickle #used to save output data to plot later...
import click
import numpy as np
import numpy.polynomial.polynomial as poly
import matplotlib.pyplot as plt
import matplotlib.font_manager
import os
from csv import reader
import sympy as sy
from scipy import signal
from pathlib import Path
import pyfiglet
import datetime
from time import sleep

def f(question):
    """Answer to the Ultimate Question of Life, the Universe, and Everything.
    tbd: Find the inverse function..."""  
    return 42


def read_toml(file):
    """This function loads the data of the vehicle."""
    vehicle = toml.load(file)
    Felge_d = vehicle['basic']['rim_dia']*0.0254
    H = vehicle['basic']['H_B'] * vehicle['basic']['B'] * 0.01
    vehicle['basic']['wheel_radius'] = (Felge_d + (2*H))/2
    vehicle['basic']['i_g'] = vehicle['basic']['n_max'] / 60 *(vehicle['basic']['wheel_radius'] * 2 * np.pi) / (vehicle['basic']['v_max']/3.6)

    return(vehicle)

def read_csv(file):
    """Read csv export from automeris.io
    Column Separator: comma [,]
    decimal separator: dot[.]"""
    with open(file, 'r') as read_obj: # read csv file as a list of lists
        csv_reader = reader(read_obj) # pass the file object to reader() to get the reader object
        list_of_rows = list(csv_reader) # Pass reader object to list() to get a list of lists
        Data = np.asarray(list_of_rows).astype(float)
        return Data

def rolling_fiction(vehicle, v, alpha):
    """This function calculates the rolling friction acording to Mitschke, Wallentowitz 2014 - Dynamik der Kraftfahzeuge p. 14.:
    
    .. math::
        
        F_{R} = (M_{w} + M_{z}) \\cdot g \\cdot cos(\\alpha) \\cdot \\left( f_{R0} + f_{R1} \\cdot \\frac{v}{100 km/h} + f_{R4} \\cdot  {\\left(\\frac{v}{100 km/h}\\right)}^4\\right)
    
    
    """
    if v == 0:
        return(0)
    else:
        return((vehicle['basic']['m_w'] + vehicle['basic']['m_z']) * vehicle['basic']['g_earth'] * np.cos(alpha /180 * np.pi) * (vehicle['basic']['f_B'] + vehicle['basic']['f_R0'] + (vehicle['basic']['f_R1'] * (v*3.6)/100) + (vehicle['basic']['f_R4'] * ((v*3.6)/100)**4)))

def aerodynamic_resistance(vehicle, v):
    """This function calculates the aerodynamic resistance acording to:
    
    .. math::
        
        F_{L} = c_{w} \\cdot A_{front} \\cdot \\frac{p_{air}}{2} \\cdot v^2
    
    
    """
    return(0.5 * vehicle['basic']['cw'] * vehicle['basic']['A_front'] * vehicle['basic']['p_air'] * v**2)

def acceleration_resistance(a, vehicle):
    """This function calculates the acceleration resistance acording to:
    
   .. math::
        
        F_{A} = a \\cdot (M_{w} \\cdot \\lambda + M_{z})
     
    
    """
    return(a * (vehicle['basic']['m_w'] * vehicle['basic']['lambda'] + vehicle['basic']['m_z']))

def slope_resistance(vehicle, alpha):
    """This function calculates the slope resistance acording to:
    
    .. math::
    
        F_{S} = sin (\\alpha) \\cdot (M_{w} + M_{z}) \\cdot g
    
    """
    return(np.sin (alpha * np.pi/180) * (vehicle['basic']['m_w'] + vehicle['basic']['m_z']) * vehicle['basic']['g_earth'])

def total_force(vehicle, a, v, alpha):
    """This function calculates the total force for a vehicle in one operation point:
    
    .. math::
    
        F_{Total} = F_{S} + F_{A} + F_{L} + F_{R}
        
    """
    return(slope_resistance(vehicle, alpha) + acceleration_resistance(a, vehicle) + aerodynamic_resistance(vehicle, v) + rolling_fiction(vehicle, v, alpha))

def dq0_to_abc(d, q, wt, z=0, delta=0):
    """Inverse Park transform"""
    a = d*np.cos(wt+delta) - q*np.sin(wt+delta) + z
    b = d*np.cos(wt-(2*np.pi/3)+delta) - q*np.sin(wt-(2*np.pi/3)+delta) + z
    c = d*np.cos(wt+(2*np.pi/3)+delta) - q*np.sin(wt+(2*np.pi/3)+delta) + z
    return a, b, c

def calc_PMSM(vehicle, T, n, verbose = False):
    """Calculation of the inverter set point for a given vehicle, tourque and speed of a permanent magnet synchronos motor (PMSM).
    ToDo: speed up -> add complete numerical calculation
    """
    erg = {}
    if T == 0:
        erg["status"] = 'Status: Basic Speed Range'
        erg["i_d"] = 0
        erg["i_q"] = 0
        erg["v_d"] = 0
        erg["v_q"] = 0
        return (erg)
     
 
    if n > vehicle['basic']['n_max']/60:
        erg["status"] = 'Error: Rotational limit reached!'
        return(erg)

    i_d = sy.symbols('i_d')
    i_dlist = (sy.solveset(sy.Eq(sy.diff(( i_d**2 + ( (2*T) / ( 3* vehicle['electric']['zp'] * (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*i_d) ) )**2) ,i_d),0), i_d, sy.Interval(-vehicle['electric']['I_max'],0)))
    if len(list(i_dlist)) == 0:
        erg["status"] = 'Error: Current limit reached!'
        return(erg)
    i_ds = list(i_dlist)[0]
    erg["i_d"] = i_ds
    i_q = (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])* i_ds ))
    erg["i_q"] = i_q
    if verbose:
        erg["fuc_T"] = sy.lambdify(i_d, (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*i_d )) )
        erg["MTPC"] = (i_ds, i_q)
    if i_ds**2 + i_q**2  > vehicle['electric']['I_max']**2:
        erg["status"] = 'Error: Current limit reached!'
        return(erg)
    omega = 2*sy.pi* n*vehicle['electric']['zp']
    v_max = vehicle['electric']['V_DC'] / np.sqrt(3) * vehicle['electric']['V_DC_res']
    v_d = (i_ds * vehicle['electric']['Rs'] - omega*vehicle['electric']['Lq'] * i_q).evalf()
    erg["v_d"] = v_d
    v_q = (i_q * vehicle['electric']['Rs'] + (vehicle['electric']['Ld'] * i_ds + vehicle['electric']['PSI'])* omega).evalf()
    erg["v_q"] = v_q
    
    if verbose:
        if T > 0: 
            erg['fuc_Imax'] = sy.lambdify(i_d, sy.sqrt(vehicle['electric']['I_max']**2 - i_d**2))
        elif T < 0:
            erg['fuc_Imax'] = sy.lambdify(i_d, -sy.sqrt(vehicle['electric']['I_max']**2 - i_d**2))

    if sy.sqrt( v_d**2 + v_q**2) > v_max:
        if T > 0:                
            func = sy.lambdify(i_d, (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*i_d ) ) - sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ))
            with np.errstate(invalid='ignore'):
                xx = np.linspace (-vehicle['electric']['I_max'],0,10000)
                if np.nanmin(func(xx)) > 0:
                    # erg["fuc_V"] = func
                    erg["fuc_V"] = sy.lambdify(i_d, sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ) )
                    erg["status"] = 'Error: Voltage limit reached!'
                    return(erg)                    

            with np.errstate(invalid='ignore'):
                A = np.diff(np.sign(func(xx)))
            A[np.isnan(A)] = 0
            test = xx[np.where(A)[0]]
            
            i_ds = test[0]
            erg["i_d"] = i_ds
            i_q = ( (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*list(test)[0] )) ) #.evalf() )
            erg["i_q"] = i_q
            erg["status"] = 'Status: Field-Weakening!'
            v_d = (i_ds * vehicle['electric']['Rs'] - omega*vehicle['electric']['Lq'] * i_q).evalf()
            erg["v_d"] = v_d
            v_q = (i_q * vehicle['electric']['Rs'] + (vehicle['electric']['Ld'] * i_ds + vehicle['electric']['PSI'])* omega).evalf()
            erg["v_q"] = v_q
            if verbose:
                erg["fuc_V"] = sy.lambdify(i_d, sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ) )
            if i_ds**2 + i_q**2  > vehicle['electric']['I_max']**2:
                erg["status"] = 'Error: Current limit under voltage limit reached!'
                return(erg)
            return(erg)
        elif T < 0:                
            func = sy.lambdify(i_d, (2*T) / ( 3*vehicle['electric']['zp']* (vehicle['electric']['PSI'] + (vehicle['electric']['Ld']-vehicle['electric']['Lq'])*i_d ) ) + sy.sqrt( (v_max**2 - (omega * vehicle['electric']['PSI'] + omega * vehicle['electric']['Ld']*i_d)**2)  / ( omega**2 * vehicle['electric']['Lq']**2) ))
            with np.errstate(invalid='ignore'):
                xx = np.linspace (-vehicle['electric']['I_max'],0,10000)