Source code for femagtools.machine

import logging
import numpy as np
import numpy.linalg as la
from .bch import Reader

import scipy.optimize as so
import scipy.interpolate as ip


logger = logging.getLogger(__name__)


[docs]def mesh(x, y): """return the combined vectors x and y""" size = np.asarray(x).size, np.asarray(y).size return (np.repeat(x, size[1]), np.tile(y, size[0]))
[docs]def K(d): """space phasor transformation matrix (Inverse Park Transformation) T-1 * dq arguments: d: rotation angle returns transformation matrix""" return np.array(( (-np.cos(d), np.sin(d)), (-np.cos(d-2*np.pi/3), np.sin(d-2*np.pi/3)), (-np.cos(d+2*np.pi/3), np.sin(d+2*np.pi/3))))
[docs]def T(d): """space phasor transformation matrix (Park Transformation) T * abc arguments: d: rotation angle returns transformation matrix""" return np.array(( (-np.cos(d), -np.cos(d-2*np.pi/3), -np.cos(d+2*np.pi/3)), (np.sin(d), np.sin(d-2*np.pi/3), np.sin(d+2*np.pi/3))))/3*2
[docs]def invpark(a, q, d): """ convert a dq vector to the abc reference frame (inverse park transformation) Args: a: rotation angle d: value in direct axis q: value in quadrature axis """ if np.isscalar(a) and np.isscalar(q) and np.isscalar(d): return np.dot(K(a), (q, d)) if np.isscalar(q) and np.isscalar(d): return np.array([K(x).dot((q, d)) for x in a]).T return np.array([K(x).dot((y, z)) for x, y, z in zip(a, d, q)]).T
[docs]def betai1(iq, id): """return beta and amplitude of dq currents""" return (np.arctan2(id, iq), la.norm((id, iq), axis=0)/np.sqrt(2.0))
[docs]def iqd(beta, i1): """return qd currents of beta and amplitude""" return np.sqrt(2.0)*i1*np.array([np.cos(beta), np.sin(beta)])
[docs]def create(bch, r1, ls, lfe=1, wdg=1): """create PmRelMachine from BCH Arguments: bch: BchReader or Erg object r1: winding resistance ls: winding leakage lfe: scale factor length wdg: scale factor number of windings """ m = 3 if isinstance(bch, Reader): p = bch.machine['p'] if bch.type.lower().find('psid-psiq-identification') >= 0: id = np.array(bch.psidq['id'])/wdg iq = np.array(bch.psidq['iq'])/wdg psid = wdg*lfe*np.array(bch.psidq['psid']) psiq = wdg*lfe*np.array(bch.psidq['psiq']) return PmRelMachinePsidq(m, p, psid, psiq, r1*lfe*wdg**2, id, iq, ls*wdg**2) if bch.type.lower().find('ld-lq-identification') >= 0: beta = bch.ldq['beta'] i1 = np.array(bch.ldq['i1'])/wdg psid = wdg*lfe*np.array(bch.ldq['psid']) psiq = wdg*lfe*np.array(bch.ldq['psiq']) return PmRelMachineLdq(m, p, psid=psid, psiq=psiq, r1=r1*lfe*wdg**2, i1=i1, beta=beta, ls=ls*wdg**22) raise ValueError("Unsupported BCH type {}".format(bch.type)) # must be ERG type: p = int(round(np.sqrt(2)*bch['M_sim'][-1][-1]/( m*bch['Psi_d'][-1][-1] * bch['i1'][-1]))) return PmRelMachineLdq(m, p, r1=r1*lfe*wdg**2, beta=bch['beta'], i1=np.array(bch['i1'])/wdg, psid=wdg*lfe*np.array(bch['Psi_d'])/np.sqrt(2), psiq=wdg*lfe*np.array(bch['Psi_q'])/np.sqrt(2), ls=ls*wdg**2)
[docs]class PmRelMachine(object): """Abstract base class for PmRelMachines Args: m: number of winding phases p: number of pole pairs r1: stator winding resistance (in Ohm) """ def __init__(self, m, p, r1, ls): self.p = p self.m = m self.r1 = r1 self.ls = ls self.io = (1, -1)
[docs] def torque_iqd(self, iq, id): "torque at q-d-current" psid, psiq = self.psi(iq, id) tq = self.m*self.p/2*(psid*iq - psiq*id) return tq
[docs] def iqd_torque(self, torque): """return minimum d-q-current for torque""" res = so.minimize(lambda iqd: la.norm(iqd), self.io, method='SLSQP', constraints=({'type': 'eq', 'fun': lambda iqd: self.torque_iqd(*iqd) - torque})) return res.x
[docs] def uqd(self, w1, iq, id): """return uq, ud of frequency w1 and d-q current""" psid, psiq = self.psi(iq, id) uqd = (self.r1*iq + w1*(self.ls*id + psid), self.r1*id - w1*(self.ls*iq + psiq)) logger.debug('beta i1 %s u1 %f', betai1(iq, id), la.norm(uqd)) return uqd
[docs] def w1_umax(self, u, iq, id): """return frequency w1 at given voltage u and id, iq current Keyword arguments: u -- the maximum voltage (RMS) iq, id -- the d-q currents""" w10 = np.sqrt(2)*u/la.norm(self.psi(iq, id)) return so.fsolve(lambda w1: la.norm(self.uqd(w1, iq, id))-u*np.sqrt(2), w10)[0]
[docs] def w1_u(self, u, iq, id): """return frequency w1 at given voltage u and id, iq current (obsolete, use w1_umax)""" return self.w1_umax(u, iq, id)
[docs] def w1max(self, u, iq, id): """return max frequency w1 at given voltage u and d-q current (obsolete, use w1_umax)""" return self.w1_umax(u, iq, id)
[docs] def w2_imax_umax(self, imax, umax): """return frequency at current and voltage""" return so.fsolve( lambda x: la.norm(self.mtpv(x, umax)[:2] - self.iqd_imax_umax(imax, x, umax)), np.sqrt(2)*umax/la.norm(self.psi(*self.io)))[0]
[docs] def beta_u(self, w1, u, i1): "beta at given frequency, voltage and current" return so.fsolve(lambda b: la.norm(self.uqd(w1, *(iqd(b, i1))))-u*np.sqrt(2), np.arctan2(self.io[1], self.io[0]))[0]
[docs] def iq_u(self, w1, u, id): "iq at given frequency, voltage and id current" return so.fsolve(lambda iq: la.norm(self.uqd(w1, iq, id))-u*np.sqrt(2), 0)[0]
[docs] def iqd_uqd(self, w1, uq, ud): "return iq, id current at given frequency, voltage" return so.fsolve(lambda iqd: np.array((uq, ud)) - self.uqd(w1, *iqd), (0, self.io[1]))
[docs] def i1_torque(self, torque, beta): "return i1 current with given torque and beta" i1, info, ier, mesg = so.fsolve( lambda i1: self.torque_iqd(*iqd(beta, i1))-torque, self.io[0], full_output=True) if ier == 1: return i1 raise ValueError("no solution found for torque {}, beta {}".format( torque, beta))
[docs] def i1_voltage(self, w1, u1, beta): "return i1 current with given w1, u1 and beta" i1, info, ier, mesg = so.fsolve( lambda i1: la.norm(self.uqd(w1, *iqd(beta, i1)))-np.sqrt(2)*u1, la.norm(self.io), full_output=True) if ier == 1: return i1 raise ValueError("no solution found for w1 {}, u1 {}, beta {}".format( w1, u1, beta))
[docs] def id_torque(self, torque, iq): "return d current with given torque and d-current" i0 = iqd(*self.io)[1] return so.fsolve(lambda id: self.torque_iqd(iq, id)-torque, i0)[0]
[docs] def iqd_torque_umax(self, torque, w1, u1max): "return d-q current and torque at stator frequency and max voltage" iq, id = self.iqd_torque(torque) # check voltage if la.norm(self.uqd(w1, iq, id)) <= u1max*np.sqrt(2): return (iq, id) # decrease psi (flux weakening mode), let i1 == i1max return so.fsolve( lambda iqd: (la.norm(self.uqd(w1, *iqd)) - u1max*np.sqrt(2), self.torque_iqd(*iqd) - torque), (iq, id))
[docs] def iqd_imax_umax(self, i1max, w1, u1max): """return d-q current at stator frequency and max voltage and max current""" beta0 = self.betarange[0] beta1 = np.sum(self.betarange)/2 u0 = la.norm(self.uqd(w1, *iqd(beta0, i1max))) u1 = la.norm(self.uqd(w1, *iqd(beta1, i1max))) du = (u0 - u1) db = (beta0 - beta1) beta0 = beta0 + db/du*(u1max*np.sqrt(2) - u0) if self.betarange[0] > beta0 or self.betarange[1] < beta0: beta0 = self.betarange[0] beta, info, ier, mesg = so.fsolve( lambda b: la.norm( self.uqd(w1, *iqd(b, i1max))) - u1max*np.sqrt(2), beta0, full_output=True) if ier == 1: return iqd(beta[0], i1max) raise ValueError( "no solution found for imax {}, w1 {}, u1max {}".format( i1max, w1, u1max))
[docs] def mtpa(self, i1): """return iq, id, torque at maximum torque of current i1""" bopt, fopt, iter, funcalls, warnflag = so.fmin( lambda x: -self.torque_iqd(*iqd(x, i1)), 0, full_output=True, disp=0) iq, id = iqd(bopt[0], i1) return [iq, id, -fopt]
[docs] def mtpv(self, w1, u1): """return d-q-current, torque for voltage and frequency with maximum torque""" res = so.minimize( lambda iqd: -self.torque_iqd(*iqd), (0, self.io[1]), method='SLSQP', constraints=( {'type': 'ineq', 'fun': lambda iqd: np.sqrt(2)*u1 - la.norm(self.uqd(w1, *iqd))})) return res.x[0], res.x[1], -res.fun # la.norm(self.uqd(w1, *(res.x))))
[docs] def characteristics(self, T, n, u1max, nsamples=36): """calculate torque speed characteristics. return dict with list values of id, iq, n, T, ud, uq, u1, i1, beta, gamma, phi, cosphi, pmech, n_type Keyword arguments: T -- the maximum torque or the list of torque values in Nm n -- the maximum speed or the list of speed values in 1/s u1max -- the maximum voltage in V rms nsamples -- (optional) number of speed samples""" r = dict(id=[], iq=[], uq=[], ud=[], u1=[], i1=[], T=[], losses=[], beta=[], gamma=[], phi=[], cosphi=[], pmech=[], n=[]) if np.isscalar(T): iq, id = self.iqd_torque(T) i1max = betai1(iq, id)[1] w1 = self.w1max(u1max, iq, id) nmax = max(w1, self.w1max(u1max, *self.iqdmin(i1max)))/2/np.pi/self.p n1 = min(w1/2/np.pi/self.p, nmax) r['n_type'] = n1 logger.info("Type speed %f n: %f nmax %f", 60*n1, 60*n, 60*nmax) try: w1 = self.w2_imax_umax(i1max, u1max) n2 = self.w2_imax_umax(i1max, u1max)/2/np.pi/self.p iqmtpv, idmtpv, tq = self.mtpv(w1, u1max) if self._inrange((iqmtpv, idmtpv)): n2 = w1/2/np.pi/self.p n3 = max(n, n2) else: n3 = min(nmax, n) n2 = n3 logger.info("n1: %f n2: %f n3 : %f ", 60*n1, 60*n2, 60*n3) except ValueError: n3 = min(nmax, n) n2 = n3 speedrange = sorted( list(set([nx for nx in [n1, n2, n3] if nx <= n3]))) n1 = speedrange[0] n3 = speedrange[-1] if n2 > n3: n2 = n3 logger.info("Speed intervals %s", [60*nx for nx in speedrange]) if len(speedrange) > 2: nsamples = nsamples - int(speedrange[1]/(n3/nsamples)) dn = (n3-speedrange[1])/nsamples else: dn = n3 / nsamples for nx in np.linspace(0, n1, int(n1/dn)): r['id'].append(id) r['iq'].append(iq) r['n'].append(nx) r['T'].append(T) if nx < n3: for nx in np.linspace(nx+dn/2, n2, int(n2)//dn): w1 = 2*np.pi*nx*self.p try: iq, id = self.iqd_imax_umax(i1max, w1, u1max) except ValueError: logger.warn("ValueError at speed %f", 60*nx) break r['id'].append(id) r['iq'].append(iq) r['n'].append(nx) r['T'].append(self.torque_iqd(iq, id)) if nx < n3: for nx in np.linspace(nx+dn/2, n3, int(n3)//dn): w1 = 2*np.pi*nx*self.p try: iq, id, tq = self.mtpv(w1, u1max) if not self._inrange((iq, id)): break except ValueError: logger.warn("ValueError at speed %f", 60*nx) break r['id'].append(id) r['iq'].append(iq) r['n'].append(nx) r['T'].append(tq) logger.info("Max speed %f", 60*nx) else: for t, nx in zip(T, n): w1 = 2*np.pi*nx*self.p iq, id = self.iqd_torque_umax(t, w1, u1max) r['id'].append(id) r['iq'].append(iq) tq = self.torque_iqd(iq, id) r['T'].append(tq) r['n'].append(nx) for nx, iq, id in zip(r['n'], r['iq'], r['id']): w1 = 2*np.pi*nx*self.p uq, ud = self.uqd(w1, iq, id) r['uq'].append(uq) r['ud'].append(ud) r['u1'].append(la.norm((ud, uq))/np.sqrt(2.0)) r['i1'].append(la.norm((id, iq))/np.sqrt(2.0)) r['beta'].append(np.arctan2(id, iq)/np.pi*180.) r['gamma'].append(np.arctan2(ud, uq)/np.pi*180.) r['phi'].append(r['beta'][-1] - r['gamma'][-1]) r['cosphi'].append(np.cos(r['phi'][-1]/180*np.pi)) for nx, tq in zip(r['n'], r['T']): r['pmech'].append((2*np.pi*nx*tq)) return r
[docs] def i1beta_characteristics(self, n_list, i1_list, beta_list, u1max): """calculate i1-beta characteristics""" r = dict(id=[], iq=[], uq=[], ud=[], u1=[], i1=[], T=[], beta=[], gamma=[], phi=[], cosphi=[], pmech=[], n=[]) for n, i1, beta in zip(n_list, i1_list, beta_list): w1 = 2*np.pi*n*self.p beta = beta/180*np.pi iq, id = iqd(beta, i1) uq, ud = self.uqd(w1, iq, id) u1 = la.norm((ud, uq))/np.sqrt(2) if u1 > u1max: logger.debug("u1 %s > %s", u1, u1max) beta = self.beta_u(w1, u1max, i1) logger.debug("beta %s", beta*180/np.pi) iq, id = iqd(beta, i1) logger.debug("beta %s id, %s iq %s", beta*180/np.pi, id, iq) uq, ud = self.uqd(w1, iq, id) u1 = la.norm((ud, uq))/np.sqrt(2) logger.debug("ud %s uq %s --> u1 %s", ud, uq, u1) tq = self.torque_iqd(iq, id) r['id'].append(id) r['iq'].append(iq) r['uq'].append(uq) r['ud'].append(ud) r['u1'].append(u1) r['i1'].append(la.norm((id, iq))/np.sqrt(2)) r['T'].append(tq) r['beta'].append(np.arctan2(id, iq)/np.pi*180.) r['gamma'].append(np.arctan2(ud, uq)/np.pi*180.) r['n'].append(n) r['phi'].append(r['beta'][-1]-r['gamma'][-1]) r['cosphi'].append(np.cos(r['phi'][-1]/180*np.pi)) r['pmech'].append(w1/self.p*r['T'][-1]) return r
def _inrange(self, iqd): i1 = np.linalg.norm(iqd)/np.sqrt(2) iqmin, idmin = self.iqdmin(i1) iqmax, idmax = self.iqdmax(i1) return iqmin <= iqd[0] <= iqmax and idmin <= iqd[1] <= idmax
[docs]class PmRelMachineLdq(PmRelMachine): """Standard set of PM machine given by i1,beta parameters: p number of pole pairs m number of phases psim flux in Vs (RMS) ld d-inductance in lq q-inductance in H r1 stator resistance ls stator leakage inductance in H beta angle i1 vs up in degrees i1 current in A (RMS) optional keyword args: psid D-Flux in Vs (RMS) psiq Q-Flux in Vs (RMS) """ def __init__(self, m, p, psim=[], ld=[], lq=[], r1=0, beta=[], i1=[], ls=0, **kwargs): super(self.__class__, self).__init__(m, p, r1, ls) self.psid = None self.betarange = (-np.pi, np.pi) self.i1range = (0, np.inf) if np.isscalar(ld): self.ld = lambda b, i: ld self.psim = lambda b, i: psim self.lq = lambda b, i: lq logger.debug("ld %s lq %s psim %s", ld, lq, psim) return if len(ld) == 1: try: self.io = iqd(min(beta)*np.pi/360, max(i1)/2) except: self.io = (1, -1) self.ld = lambda b, i: ld[0] self.psim = lambda b, i: psim[0] self.lq = lambda b, i: lq[0] logger.debug("ld %s lq %s psim %s", ld, lq, psim) return beta = np.asarray(beta)/180.0*np.pi self.io = iqd(np.min(beta)/2, np.max(i1)/2) if 'psid' in kwargs: kx = ky = 3 if len(i1) < 4: ky = len(i1)-1 if len(beta) < 4: kx = len(beta)-1 self.betarange = min(beta), max(beta) self.i1range = (0, np.max(i1)) psid = np.sqrt(2)*np.asarray(kwargs['psid']) psiq = np.sqrt(2)*np.asarray(kwargs['psiq']) self.psid = lambda x, y: ip.RectBivariateSpline( beta, i1, psid, kx=kx, ky=ky).ev(x, y) self.psiq = lambda x, y: ip.RectBivariateSpline( beta, i1, psiq, kx=kx, ky=ky).ev(x, y) return if len(i1) < 4 or len(beta) < 4: if len(i1) == len(beta): self.ld = lambda x, y: ip.interp2d(beta, i1, ld.T)(x, y) self.psim = lambda x, y: ip.interp2d(beta, i1, psim.T)(x, y) self.lq = lambda x, y: ip.interp2d(beta, i1, lq.T)(x, y) logger.debug("interp2d beta %s i1 %s", beta, i1) return elif len(i1) == 1: self.ld = lambda x, y: ip.InterpolatedUnivariateSpline( beta, ld, k=1)(x) self.psim = lambda x, y: ip.InterpolatedUnivariateSpline( beta, psim, k=1)(x) self.lq = lambda x, y: ip.InterpolatedUnivariateSpline( beta, lq, k=1)(x) logger.debug("interpolatedunivariatespline beta %s", beta) return if len(beta) == 1: self.ld = lambda x, y: ip.InterpolatedUnivariateSpline( i1, ld, k=1)(y) self.psim = lambda x, y: ip.InterpolatedUnivariateSpline( i1, ld, k=1)(y) self.lq = lambda x, y: ip.InterpolatedUnivariateSpline( i1, lq, k=1)(y) logger.debug("interpolatedunivariatespline i1 %s", i1) return raise ValueError("unsupported array size {}x{}".format( len(beta), len(i1))) self.betarange = min(beta), max(beta) self.i1range = (0, np.max(i1)) self.ld = lambda x, y: ip.RectBivariateSpline( beta, i1, np.asarray(ld)).ev(x, y) self.psim = lambda x, y: ip.RectBivariateSpline( beta, i1, np.asarray(psim)).ev(x, y) self.lq = lambda x, y: ip.RectBivariateSpline( beta, i1, np.asarray(lq)).ev(x, y) logger.debug("rectbivariatespline beta %s i1 %s", beta, i1)
[docs] def psi(self, iq, id): """return psid, psiq of currents iq, id""" beta, i1 = betai1(np.asarray(iq), np.asarray(id)) logger.debug('beta %f (%f, %f) i1 %f %f', beta, self.betarange[0], self.betarange[1], i1, self.i1range[1]) if (self.betarange[0] <= beta <= self.betarange[1] and i1 <= 1.01*self.i1range[1]): if self.psid: return (self.psid(beta, i1), self.psiq(beta, i1)) psid = self.ld(beta, i1)*id + np.sqrt(2)*self.psim(beta, i1) psiq = self.lq(beta, i1)*iq return (psid, psiq) return (np.nan, np.nan)
[docs] def psi0(self, iq, id): """return psid, psiq of currents iq, id""" beta, i1 = betai1(np.asarray(iq), np.asarray(id)) if self.psid: return (self.psid(beta, i1), self.psiq(beta, i1)) psid = self.ld(beta, i1)*id + np.sqrt(2)*self.psim(beta, i1) psiq = self.lq(beta, i1)*iq return (psid, psiq)
[docs] def iqdmin(self, i1): """max iq, min id for given current""" return iqd(self.betarange[0], i1)
[docs] def iqdmax(self, i1): """max iq, min id for given current""" return iqd(self.betarange[1], i1)
[docs]class PmRelMachinePsidq(PmRelMachine): """Standard set of PM machine parameters: p number of pole pairs m number of phases psid d-flux (Vs Peak) psiq q-flux (Vs Peak) r1 stator resistance (Ohm) r1 stator leakage inductance (H) id q current (A, Peak) iq q current (A, Peak) """ def __init__(self, m, p, psid, psiq, r1, id, iq, ls=0): super(self.__class__, self).__init__(m, p, r1, ls) if isinstance(psid, (float, int)): self._psid = lambda id, iq: np.array([[psid]]) self._psiq = lambda id, iq: np.array([[psiq]]) return psid = np.asarray(psid) psiq = np.asarray(psiq) id = np.asarray(id) iq = np.asarray(iq) self.idrange = (min(id), max(id)) self.iqrange = (min(iq), max(iq)) if np.mean(self.idrange) > 0: self.betarange = (0, np.pi/2) else: self.betarange = (-np.pi/2, 0) self.io = np.max(iq)/2, np.min(id)/2 if np.any(psid.shape < (4, 4)): if psid.shape[0] > 1 and psid.shape[1] > 1: self._psid = ip.interp2d(iq, id, psid.T) self._psiq = ip.interp2d(iq, id, psiq.T) return if len(id) == 1 or psid.shape[1] == 1: self._psid = lambda x, y: ip.InterpolatedUnivariateSpline( iq, psid)(x) self._psiq = lambda x, y: ip.InterpolatedUnivariateSpline( iq, psiq)(x) return if len(iq) == 1 or psid.shape[0] == 1: self._psid = lambda x, y: ip.InterpolatedUnivariateSpline( id, psid)(y) self._psiq = lambda x, y: ip.InterpolatedUnivariateSpline( id, psiq)(y) return raise ValueError("unsupported array size {}x{}".format( len(psid.shape[0]), psid.shape[1])) self._psid = lambda x, y: ip.RectBivariateSpline( iq, id, psid).ev(x, y) self._psiq = lambda x, y: ip.RectBivariateSpline( iq, id, psiq).ev(x, y)
[docs] def psi(self, iq, id): return (self._psid(iq, id), self._psiq(iq, id))
[docs] def iqdmin(self, i1): """max iq, min id for given current""" if np.min(self.idrange) < 0 and np.max(self.idrange) <= 0: idmin = -np.sqrt(2)*i1 else: idmin = 0 if np.min(self.idrange) <= idmin: iqmin = 0 if np.min(self.iqrange) <= iqmin: return (iqmin, idmin) return np.min(self.iqrange), idmin i1max = np.sqrt(2)*i1 beta = np.arccos(self.iqrange[0]/i1max) iqmin = np.sqrt(2)*i1*np.sin(beta) if np.min(self.iqrange) <= iqmin: return (iqmin, idmin) return np.min(self.iqrange), np.min(self.idrange)
[docs] def iqdmax(self, i1): """max iq, max id for given current""" iqmax = np.sqrt(2)*i1 if iqmax <= np.max(self.iqrange): if np.min(self.idrange) < 0 and np.max(self.idrange) <= 0: idmax = 0 else: idmax = np.sqrt(2)*i1 if idmax <= np.max(self.idrange): return (iqmax, idmax) return (iqmax, np.max(self.idrange)) beta = np.arccos(self.iqrange[1]/iqmax) iqmax = np.sqrt(2)*i1*np.cos(beta) idmax = np.sqrt(2)*i1*np.sin(beta) if idmax <= np.max(self.idrange): return (iqmax, idmax) return iqmax, np.max(self.idrange)