#!/usr/bin/python3 from math import * import numpy as np from scipy.fft import fft, ifft, fftfreq, next_fast_len from scipy.optimize import brentq as findRoot import subprocess from time import sleep ######################################################################################################### class gnuPlotter: def __init__(self): self.gpp={} def open_gpp(self, gppnamestr=None): if (gppnamestr is None) or (gppnamestr not in self.gpp.keys()): gppnew=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE) #gppnew=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE, stderr=subprocess.PIPE) gppnew.stdin.write("set term qt size 1900,950; \n") #gppnew.stdin.write("set term x11 size 1900,950; \n") if gppnamestr is None: return gppnew if gppnamestr not in self.gpp.keys(): self.gpp[gppnamestr]=gppnew return self.gpp[gppnamestr] #def gpp_Raw(self, ldat, rawstr, gppobj=None): def gpp_Raw(self, ldat, rawstr, gppstr=None): gppCurr=self.open_gpp(gppstr) #write ldat to datablock gppCurr.stdin.write("$idat << EOD \n") for irow in ldat: gppCurr.stdin.write(("{} "*len(irow)+" \n").format(*irow)) gppCurr.stdin.write("EOD \n") #send raw command; if datafile/block is passed to script, use '$idat' gppCurr.stdin.write(rawstr) #def gpp_Plot(self, ldat, strdat, lstrlbl, gppobj=None, lineWidth=2, toFile=''): def gpp_Plot(self, ldat, strdat, lstrlbl, gppstr=None, lineWidth=2, toFile='', config=''): gppCurr=self.open_gpp(gppstr) gppCurr.stdin.write(config) lcol=[r.split("[")[0].strip() for r in strdat.split(",")] def get_plotStr(strlbl): llbl=[l.strip().split(":") for l in strlbl.split(",")] nlbl=len(llbl) lidx=[ [ [lcol.index(s) for s in lcol if s.endswith(lbl[ii])] for lbl in llbl] for ii in [0,1] ] lidxpairs=[(lidx[0][i][0],lidx[1][i][0]) for i in range(0,nlbl)] listr=["{}:{}".format(i[0]+1,i[1]+1) for i in lidxpairs] plotStr="plot " #litemStr=["$idat u {} w l lw 2 title '{}'".format(listr[i], llbl[i][1]) for i in range(0,len(listr))] litemStr=["$idat u {} w l lw {} title '{}'".format(listr[i], lineWidth, llbl[i][1]) for i in range(0,len(listr))] return plotStr+", ".join(litemStr)+" \n" lplotStr=[get_plotStr(lbl) for lbl in lstrlbl] if toFile.endswith(".pdf"): gppCurr.stdin.write("set term pdfcairo enhanced font ',6' \n") #print("set output '{}' \n".format(toFile)) gppCurr.stdin.write("set output '{}' \n".format(toFile)) #write ldat to datablock gppCurr.stdin.write("$idat << EOD \n") for irow in ldat: gppCurr.stdin.write(("{} "*len(irow)+" \n").format(*irow)) gppCurr.stdin.write("EOD \n") #plot or multiplot if len(lstrlbl)==1: gppCurr.stdin.write(lplotStr[0]) elif len(lstrlbl)==2: gppCurr.stdin.write("set multiplot \n") gppCurr.stdin.write("set origin 0,0 \n set size 1,0.5 \n") gppCurr.stdin.write(lplotStr[0]) gppCurr.stdin.write("set origin 0,0.5 \n set size 1,0.5 \n") gppCurr.stdin.write(lplotStr[1]) gppCurr.stdin.write("unset multiplot \n") #gppCurr.stdin.write(lplotStr[1]+" pause -1 \n") elif len(lstrlbl)==3: gppCurr.stdin.write("set multiplot \n") gppCurr.stdin.write("set origin 0,0 \n set size 1,0.35 \n") gppCurr.stdin.write(lplotStr[0]) gppCurr.stdin.write("set origin 0,0.333 \n set size 1,0.35 \n") gppCurr.stdin.write(lplotStr[1]) gppCurr.stdin.write("set origin 0,0.666 \n set size 1,0.35 \n") gppCurr.stdin.write(lplotStr[2]) gppCurr.stdin.write("unset multiplot \n") #gppCurr.stdin.write(lplotStr[1]+" pause -1 \n") #if toFile is not None: if toFile != '': gppCurr.stdin.write("set output \n") return gppCurr ######################################################################################################### class fftIntegrator: def __init__(self, nptspp, fperiod, npad, nprev, ztype2fn, ztype2fn_dc, resultfn, resultbounds): self.nPtsPP=nptspp self.fPeriod=fperiod self.nPtsPS=round(self.nPtsPP*self.fPeriod) self.tauLen=1.0 self.tauScale=1.0/self.fPeriod self.nPad=npad self.nPrev=nprev self.nPerTot=self.nPad+self.nPrev+self.nPad self.nPtsTot=self.nPtsPP*self.nPerTot self.zType2Fn=ztype2fn self.zType2Fn_DC=ztype2fn_dc self.resultFn=resultfn self.resultBounds=resultbounds self.tArr=None self.freqArr=None self.prevprevArr=None self.prevArr=None self.padArr=None self.gpp={} def fftIntegrate(self, arrcurr, flagSavePrev, plotit=False): if not "fftInt" in self.gpp.keys(): self.gpp["fftInt"]=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE, stderr=subprocess.PIPE) self.gpp["fftInt"].stdin.write("set term x11 size 1900,950; \n") if self.freqArr is None: self.freqArr=fftfreq(self.nPtsTot, d=self.tauScale/(self.nPtsPP)) if self.tArr is None: tpad_pre = np.linspace(-(self.nPrev+self.nPad)*self.tauLen, -self.nPrev*self.tauLen, round(self.tauScale*self.nPad*self.nPtsPS), endpoint=False) tprev = np.linspace(-self.nPrev*self.tauLen, 0.0, round(self.tauScale*self.nPrev*self.nPtsPS), endpoint=False) tpad_post = np.linspace(0.0, self.nPad*self.tauLen, round(self.tauScale*self.nPad*self.nPtsPS), endpoint=False) #first period starting t=0 overwritten by current period to integrate self.tArr = np.hstack((tpad_pre, tprev, tpad_post)) if self.prevArr is None: self.prevArr=np.zeros(round(self.tauScale*self.nPrev*self.nPtsPS)) self.padArr=np.zeros(round(self.tauScale*self.nPad*self.nPtsPS)) idx0=-round(self.tauScale*self.nPad*self.nPtsPS) idxf=-round(self.tauScale*(self.nPad-1)*self.nPtsPS) idxshift_beg=-round(self.tauScale*(self.nPad+self.nPrev-1)*self.nPtsPS) #one period past start of previous history idxshift_end=-round(self.tauScale*(self.nPad-1)*self.nPtsPS) #to end of current solution, which=idxf if self.tArr[idx0]!=0.0 or self.tArr[idxf]!=1.0: print("fftIntegrate: ERROR: check tcurr boundaries, should be 0.0, 1.0 by definition:", idx0, self.tArr[idx0], idxf, self.tArr[idxf]) quit() ft=np.hstack((self.padArr, self.prevArr, arrcurr, self.padArr[idxf:])) result = self.fftIntegrate_Type2(ft, idx0) rr=result[idx0:idxf] vaoclip=np.less(self.resultFn()+result.real, self.resultBounds[1])*np.greater(self.resultFn()+result.real, self.resultBounds[0]) #ft= ft*vaoclip #not the most accurate method, but seems best ft= ft-0.9*ft*np.logical_not(vaoclip) result=np.clip(result.real, a_min=self.resultBounds[0]-self.resultFn(), a_max=self.resultBounds[1]-self.resultFn()) if plotit>1: self.gpp["fftInt"].stdin.write("$idat << EOD \n") for i in range(0, len(ft)): ll=[self.tArr[i], ft[i], result[i].real, int(vaoclip[i])] self.gpp["fftInt"].stdin.write(("{} "*len(ll)+" \n").format(*ll)) self.gpp["fftInt"].stdin.write("EOD \n") self.gpp["fftInt"].stdin.write("set multiplot \n") self.gpp["fftInt"].stdin.write("set origin 0,0 \n set size 1,0.5 \n") #self.gpp["fftInt"].stdin.write("plot $idat u 1:3 w l lw 2 title 'result' \n") self.gpp["fftInt"].stdin.write("plot $idat u 1:3 w l lw 2 title 'result', $idat u 1:4 w l lw 2 title 'vaoclip' \n") self.gpp["fftInt"].stdin.write("set origin 0,0.5 \n set size 1,0.5 \n") self.gpp["fftInt"].stdin.write("plot $idat u 1:2 w l lw 2 title 'ft' \n") self.gpp["fftInt"].stdin.write("unset multiplot \n") if flagSavePrev: self.prevprevArr=self.prevArr self.prevArr=ft[idxshift_beg:idxshift_end] return rr def fftIntegrate_Type2(self, ft, idx0): fw = fft(ft, norm="forward") #print("freq[1]={} freq[{}]={} self.tauScale={} t0=t[{}]={} sum(fw)={:1.6e}".format(self.freqArr[1], (self.nPtsTot/2)-1, self.freqArr[int(self.nPtsTot/2)-1:int(self.nPtsTot/2)+1], self.tauScale, idx0, self.tArr[idx0], sum(fw))) fwlp = np.hstack((0, self.zType2Fn(fw[1:], self.freqArr[1:]))) ftlp = ifft(fwlp, norm="forward") #print("DC components: fw[0]={:1.6e}, ftlp[0]={:1.6e}, self.tauScale={:1.6e}, self.zType2Fn_DC().real={:1.6e}".format(fw[0], ftlp[0], self.tauScale, self.zType2Fn_DC().real)) ftlp += (self.zType2Fn_DC().real*fw[0]*self.tauScale*(self.tArr-self.tArr[idx0])) - ftlp[idx0] return ftlp ######################################################################################################### # UCC28070 ######################################################################################################### class UCC28070: def __init__(self, plotpdf=False): ####################### # DESIGN PARAMETERS ####################### # Voltage Source AC V RMS self.Vac_rms=240 self.Fac=60 # Voltage Output DC V self.Vout=600 self.Vovp=self.Vout*1.06 #by definition self.Vdrop=20 #constraints on Cout; max Vdrop in time thalf; start with Vripple instead? # Power Total W RMS self.Ptot=1600 # Frequency Switching PWM self.Fpwm=60*3334 #use int multiple of Fac self.Fsync=2*self.Fpwm self.Dmax=0.98 ####################### # INTERNAL PARAMETERS ####################### # Regulation Voltages self.Vreg_vsense=3.0 #not adjustable self.Vreg_vinac=3.0 #Disable Synthesis --> adjustable: optimal ~2.7774 = sqrt(Kvff_max) #self.Vreg_vinac=2.7774 #Disable Synthesis --> adjustable: optimal ~2.7774 = sqrt(Kvff_max) self.DVvao=3.6 #Vvao range, 5-1=4.0 max, leave room for power loss, vac fluctuations #self.Vvao_bounds=[0.0, 5.0] self.Vreg_currsense=3.2 #max cm volts input 3.6V, leave room for peaks self.Vreg_currpeak=4.0 #must have room for ILxPeak, up to 4V self.DVramp=4.0 #not adjustable self.Vref=6.0 # Voltage Amplifier Transconductance VAO self.Gmv=70.0e-6 #uS # Current Amplifier Transconductance CAOx self.Gmc=100.0e-6 #uS # Current Multiplier IMO self.Kmult=17.0e-6 #uA ####################### # COMPONENTS ####################### # Inductors self.Lx=540.0e-6 # VSENSE divider self.RAvs=3.0e6; self.RBvs=self.get_RB(self.RAvs, self.Vreg_vsense, self.Vout) #0.2mA --> only 120mW self.Cvs=1.0e-9 # VINAC divider (div ratio should match vsense ratio, but with NO synthesis and Kvff error, we can in fact rescale vinac a bit to cancel error, best match actual power output # ??? NOT full 3V range, most useful for 120VAC applications!! more accurately detect zero crossing and better signal/noise!! hooray for no current synthesis) # BUT can I use full scale and simply rescale Rimo to compensate for rescaled kvff? I think yes... self.RAac=3.0e6; self.RBac=self.get_RB(self.RAac, self.Vreg_vinac, self.Vout) #self.RAac=3.0e6; self.RBac=self.get_RB(self.RAac, self.Vreg_vinac, self.Vac_rms*sqrt(2.0)) self.Cac=1.0e-9 # Dither / Sync self.Rrdm=0 #GND:no sync no dither; >5V:sync to pulse on RDM; self.Ccdr=0 #dither off: CDR to VREF; # IMO current mult output; Rimo=((1/2)*Iin_pk*Rsequiv)/Imo_max; Imo_max=17uA*(Vac_rms_lopk)*(Vvao-1)/kvff; kvff=0.398, Vvao=5V #self.Rimo=17500 #using 3V/171uA self.Rimo=self.get_Rimo() self.Rsequiv=self.get_Rsequiv() # Current Synthesis downslope; VRsynth = 3.0 pm0.09 V; TI info USELESS. propto L*(RB/(RA+RB))/tdrop, tdrop=1/(Rsequiv*0.1nF), Rsequiv=3V/Ipkavg self.Enable_Synthesis=0 #disable synthesis if full-time current sensing used self.Rsyn=self.get_Rsyn() # Peak Current limit; Vpklimit = 3.3 pm0.03 V at CSx rising self.RApk=6.8e3; self.RBpk=self.get_RB(self.RApk, self.Vreg_currpeak, self.Vref) #choose to draw <<1mA from Vref=6V #self.RBpk=6.6e3 #3.5 when div app to Vref=6V # Oscillator Rate; Rrt(kOhm) = 15000 / fsync(kHz); seems internal C=66.6666pF #Duty Cycle Max; Rdmx = (Rrt)*(2*Dmax-1-Dsync); Dsync=tsync/(1/fsync); tsync>=200nS; internal osc -->Dsync=0 self.PWsyncExt=250.0e-9 #external sync pulse-width (not period); minimum from datasheet -> 200.0e-9 self.Rrt, self.Rdmx = self.get_Rrt_Rdmx() # Vout filter self.Cout=self.get_Cout(self.Vdrop) # VAO compensation type 2 self.fracDVvao=0.02 #fraction of DVvao equivalent to Vripple self.Cpv=self.get_Cpv() self.Fpv=self.get_Fpv() self.Fzv=(1/6)*self.get_Fpv() #set zero freq to fraction of pole freq self.Rzv, self.Czv = self.get_Rzv_Czv() # CAOx compensation type 2 self.Fzc=(1/12)*self.Fpwm self.Gpsc=abs(self.get_Gpsc(self.Fzc)) #Gain: power stage current (magnitude) self.Rzc, self.Czc = self.get_Rzc_Czc() self.Fpc=(1/3)*self.Fpwm self.Cpc=self.get_Cpc() # Vref,Vcc filters self.Cvref=1.0e-6 #minimum 22nF self.Cvcc=10.0e-6 #minumum 0.1uF # Soft Start; recommend Css >= Czv to limit overshoot; try optimal finite time instead self.Css=self.get_Css(0.5) #300ms too fast? ####################### # STATE VARIABLES ####################### self.Clock=0 self.Rload=225.0 #ideal pure resistive load self.EnablePh1=True #False-> disable second phase (Ph1) self.SkipHalfMod=1 #1->no skip n->skip n self.SkipHalf=False #dynamic, computed in VLoop self.gppPlotter=gnuPlotter() self.PlotPDF=plotpdf self.PlotPrefix="pdfplots/" self.VL0 = self.VoltageLoop(self) self.CL0 = self.CurrentLoop(self) self.PM0 = self.PowerMeter() ######################################################################################################### # OSCILLATOR : Rate, Duty Max def get_DsyncExt(self, tsync, tpropdly=100.0e-9): return self.PWsyncExt*self.Fsync #no prop dly #return (self.PWsyncExt-tpropdly)*self.Fsync #with prop dly def get_Rrt_Rdmx(self): rrt = (15.0e9 / self.Fsync) rdmx = rrt * (2.0*self.Dmax - 1.0 - self.get_DsyncExt(self.PWsyncExt)) return rrt, rdmx def get_RB(self, ra, numer, denom): kfrac=numer/denom return (kfrac*ra)/(1.0-kfrac) def get_Css(self, tss): return tss*(10.0e-6/2.25) #10uA charging cap to 2.25V in ~tss sec ######################################################################################################### # VSENSE, VAO : Voltage Sense, Voltage Amplifier def get_KVsense(self): return self.RBvs/(self.RAvs+self.RBvs) def get_Iimo(self, vinac, vvao): return (self.Kmult)*vinac*(vvao-1.0)/self.get_Kvff() def get_Vimo(self, iimo): return iimo*self.Rimo*(1.0+float(self.EnablePh1==False))*(1.0+float(self.SkipHalf==True)) def get_Cout(self, vdrop): return float( ((1.0/(2.0*self.Fac))*self.Ptot)/((1/2)*((self.Vout)**2 - (self.Vout-vdrop)**2)) ) def get_Vripple(self): return self.get_ICout(sqrt(2)*self.Vac_rms, 2*self.Ptot)*(1/(2*pi*2*self.Fac*self.Cout)) # VAO : Voltage Amplifier Compensation def get_Cpv(self): ztmp=(self.fracDVvao*self.DVvao)/(self.get_Vripple()*self.get_KVsense()*self.Gmv) #print("ztmp",ztmp) return 1.0/(2*pi*(2*self.Fac)*ztmp) def get_Fpv(self): return (1/(2*pi))*sqrt( (self.get_KVsense()*self.Gmv*self.Ptot)/(self.DVvao*self.Vout*self.Cout*self.get_Cpv()) ) def get_Rzv_Czv(self): rzv=1/(2*pi*self.get_Cpv()*self.get_Fpv()) czv=1/(2*pi*rzv*self.Fzv) return rzv, czv # Gcv : Voltage Compensation Network Gain def get_Gcv(self, f): return self.get_KVsense() * self.Gmv * self.get_Zcomp_type2(1.0, f, self.Rzv, self.Czv, self.Cpv) # Gpsv : Voltage Loop Power Stage Gain def get_Gpsv(self, f): return ((self.Ptot/self.Vout)/(1j*2*pi*f*self.Cout))/self.DVvao # GTOTv : TOTAL VOLTAGE LOOP GAIN def get_GTOTv(self, f): return self.get_Gcv(f) * self.get_Gpsv(f) ######################################################################################################### # IMO : Current Multiplier Output def get_IMO_max(self): #if Vac increases or decreases, what is possible change in current? within kvff level? return (self.Kmult)*(self.get_Vvinac(sqrt(2)*self.Vac_rms))*(self.DVvao)/self.get_Kvff() def get_Rimo(self): return float(self.Vreg_currsense/self.get_IMO_max()) def get_KVvinac(self): return self.RBac/(self.RAac+self.RBac) def get_Vvinac(self, vinac): return vinac*self.get_KVvinac() def get_Vinac_Correction_Factor(self): vpk=self.get_Vvinac(sqrt(2)*self.Vac_rms) kvpk=self.get_Kvff() result=[vpk, kvpk, (vpk**2)/2.0, ((vpk**2)/2.0)/kvpk, sqrt(2.0*kvpk), sqrt(2.0*kvpk)/vpk] if self.PlotPDF: #print("result:",result) def fk(k): return 3.0*k/1000.0 self.kvffDat = [ [ fk(k), self._get_Kvff(fk(k)), fk(k)**2/2.0, (fk(k)**2/2.0)/self._get_Kvff(fk(k)), vpk, kvpk, (vpk**2/2.0)/(kvpk), 1.0] for k in range(1,1000) ] self.kvff_str= "Vvinac, Kvff, a*Pvinac, Pratio, VPeak, KVPeak, KVPeakRatio, Perfect" self.kvff_lblstr=["Vvinac:Kvff, Vvinac:a*Pvinac, VPeak:Kvff","Vvinac:Pratio, Vvinac:KVPeakRatio, VPeak:Pratio, Vvinac:Perfect"] self.gppPlotter.gpp_Plot(self.kvffDat, self.kvff_str, self.kvff_lblstr, toFile=self.get_FileName("kvff_ratios.pdf")) #quit() return result #Kvff discrete levels cause error in IMO current command, fix this def get_Kvff(self): # NOTE this will be irrelevant when current synthesis is disabled thanks to full-time inductor current sensing via hall effect sensors # which has huge advantages for low VAC applications, we can scale Vinac divider for full accuracy zero detect and best signal/noise return self._get_Kvff(self.get_Vvinac(sqrt(2)*self.Vac_rms)) def _get_Kvff(self, v, dvin_sign=1): #klst=[[0,0], [1.0, 0.398], [1.2, 0.600 ], [1.4, 0.839], [1.65, 1.156], [1.95, 1.604], [2.25, 2.199], [2.6, 2.922], [3.0*33, 3.857]] return np.piecewise(v, \ [v<1.0, (v>=1.0)&(v<1.2), (v>=1.2)&(v<1.4), (v>=1.4)&(v<1.65), (v>=1.65)&(v<1.95), (v>=1.95)&(v<2.25), (v>=2.25)&(v<2.6), v>=2.6], \ [0.398, 0.600, 0.839, 1.156, 1.604, 2.199, 2.922, 3.857] ) # CAOx : Current Synthesizer, Amplifiers def get_Rsequiv(self): return self.Vreg_currsense/self.get_ILxTop(self.Vac_rms, self.Ptot/2) def get_Rsyn(self): return float(self.Enable_Synthesis)*self.Lx*self.get_KVsense()/((0.1e-9)*self.Rsequiv) #appears 100pF cap inside def get_Rzc_Czc(self): rzc=1.0/(self.Gmc*abs(self.get_Gpsc(self.Fzc))) czc=1.0/(2*pi*self.Fzc*rzc) return rzc, czc def get_Cpc(self): return 1/(2*pi*self.Fpc*self.Rzc) # Gcc : Current Compensation Network Gain def get_Gcc(self, f): return self.Gmc * self.get_Zcomp_type2(1.0, f, self.Rzc, self.Czc, self.Cpc) # CAOx : Current Loop Power Stage Gain def get_Gpsc(self, f): return (self.Vout*self.Rsequiv)/(1j*2*pi*f*self.Lx*self.DVramp) # GTOTc : TOTAL CURRENT LOOP GAIN def get_GTOTc(self, f): return self.get_Gcc(f) * self.get_Gpsc(f) # AC voltage, current def get_Vinac_Sin(self, theta): # Pi Periodic, PosDef for real theta return sqrt(2)*self.Vac_rms*np.fabs(np.sin(theta)) def get_Iin_Sin(self, theta, prms): # Pi Periodic, PosDef for real theta return sqrt(2)*(prms/self.Vac_rms)*np.fabs(np.sin(theta)) # Idealized approximate quantities def get_Duty(self, vin): return (self.Vout-vin)/self.Vout def get_Iin(self, vin, pout): return pout/vin def get_ILx(self, vin, poutLx): return self.get_Iin(vin, poutLx) def get_ILxTop(self, vin, poutLx): return self.get_ILx(sqrt(2)*vin, 2*poutLx) def get_ILxPeak(self, fracsign=1): return (fracsign*self.get_dILxTop(self.Vac_rms)/2.0) + self.get_ILxTop(self.Vac_rms, self.Ptot/2) def get_dILx(self, vin, freq): return vin*(self.get_Duty(vin)/freq)/(self.Lx) def get_dILxTop(self, vin): return self.get_dILx(sqrt(2)*vin, self.Fpwm) def get_LxMin_dILx(self, vin, poutLx, freq, fracilx): return sqrt(2)*poutLx / ( (self.get_ILxTop(vin, poutLx)**2)*fracilx*freq ) def get_LxMin_CCM(self, vin, poutLx, freq): return (vin**2)/(2*poutLx*freq) def get_PoutMin_CCM(self, vin, freq): return 2*(vin**2)/(2*self.Lx*freq) def get_ICout(self, vin, ptot): return 2*self.get_ILx(vin, ptot/2)*(1.0-self.get_Duty(vin)) def _get_ILx_Vals(self, n, pout): nmax=(self.Fpwm/self.Fac)/2.0 ph_on, ph_mid = pi*((n)/nmax), pi*((n+0.5)/nmax) v, i = self.get_Vinac_Sin(ph_mid), self.get_Iin_Sin(ph_mid, pout) d = self.get_Duty(v) ph_off, ph_end = pi*((n+d)/nmax), pi*((n+1)/nmax) dt_on, dt_off = d/self.Fpwm, (1-d)/self.Fpwm iL, diL = i, (v*dt_on)/(2*self.Lx) return ph_on, ph_off, ph_end, v, i, d, dt_on, dt_off, iL-diL, iL+diL def get_FileName(self, fstr): if self.PlotPDF: return self.PlotPrefix+fstr else: return '' def plotBode_VoltageLoop(self): ftmp=[p/100.0 for p in range(0, 300, 1)]; flo=1 gptmp=[ self.get_Z_lst(flo*pow(10,f), self.get_GTOTv(flo*pow(10,f))) for f in ftmp] self.gppPlotter.gpp_Raw(gptmp, "call 'gpBode.gp' '$idat' '" + self.get_FileName("bodePlot_VoltageLoop.pdf") + "' \n") def plotBode_CurrentLoop(self): ftmp=[p/100.0 for p in range(301, 600, 1)] gptmp=[ self.get_Z_lst(pow(10,f), self.get_GTOTc(pow(10,f))) for f in ftmp] self.gppPlotter.gpp_Raw(gptmp, "call 'gpBode.gp' '$idat' '" + self.get_FileName("bodePlot_CurrentLoop.pdf") + "' \n") ######################################################################################################### class PowerMeter: def __init__(self): self.dictE={} self.dictT={} def initLbl(self, lbl): self.dictE[lbl]=0.0 self.dictT[lbl]=0.0 def add_E_T(self, lbl, e, t): if lbl not in self.dictE.keys(): self.initLbl(lbl) self.dictE[lbl] += e self.dictT[lbl] += t def Report(self): print("PowerMeter: ") for k in self.dictE.keys(): print(k, self.dictE[k]/self.dictT[k], self.dictE[k], self.dictT[k]) ######################################################################################################### class VoltageLoop: def __init__(self, parent): self.parent=parent self.Vvao_bounds=[0.0, 5.0] self.nptsIvao=20 self.fcycIvao=2*self.parent.Fac self.thalf=(1.0/self.fcycIvao) self.dt=self.thalf*(1.0/self.nptsIvao) self.Tau = np.linspace(0, np.pi, self.nptsIvao, endpoint=False) self.Vinac=self.parent.get_Vinac_Sin(self.Tau) self.gppPlotter=gnuPlotter() self.fftIntegrator_Ivao = fftIntegrator(self.nptsIvao, self.fcycIvao, 20, 120, self.get_ZType2_VAO, self.get_ZType2_VAO_DC, self.get_Vvao, self.Vvao_bounds) self.Vvao=0 self.QCss=0 self.cntThalf=0 self.Clock=0 def addClockThalf(self): self.cntThalf+=1 self.Clock+=1.0/(self.parent.Fac*2.0) def get_Vss(self): return self.QCss/self.parent.Css def get_Vsense(self, vout): return vout*self.parent.get_KVsense() def get_VsenseRef(self): return min(3.0, self.get_Vss()) def get_Ivao(self, v): #v==Vsense; note -sign since interest in current coming out of vao isat = 30.0e-6 iscp = 100.0e-6*float(self.get_Vss()>4.0) vlo, vmid, vhi = -0.42857, -0.2, 0.42857 dv = v-self.get_VsenseRef() rnew = -np.piecewise(dv, \ [(dv=vlo)&(dv=vmid)&(dv=vhi)], \ [-isat-iscp, lambda dum: dum*self.parent.Gmv-iscp, lambda dum: dum*self.parent.Gmv, isat] \ ) return rnew def get_Vvao(self): return self.Vvao def get_ZType2_VAO(self, bk, f): return self.parent.get_Zcomp_type2(bk, f, self.parent.Rzv, self.parent.Czv, self.parent.Cpv) def get_ZType2_VAO_DC(self): #hack that returns approx integration constant for DC component: 1/(cp+cz) return self.parent.get_Zcomp_type2(1, 1/(2*pi*1j), 0, self.parent.Czv, self.parent.Cpv) #----------------------------------------------------------------------------------------------------------------------------------- def stepTHALF(self, n, useCurrentLoop, plotIterDat=True, plotThalfDat=False, nskip=0): if(self.Vvao<=1.0): self.Vvao=1.0001 lvIterDat=[] self.lvIterDat_str= "self.Clock, self.parent.CL0.QCOUT.get_VCout(), 100*self.get_Vss(), deltaVao, IvaoAvg, 100*self.Vvao, self.parent.Rload, 0.5*ploadavg, PinacAvg, PcoutAvg, IcoutAvg" self.lvIterDat_lblstr=["Clock:VCout(), Clock:ploadavg, Clock:PinacAvg, Clock:PcoutAvg, Clock:Vvao, Clock:Vss()"] self.lvThalfDat_str= "Tphase, nphon*Iinac, Vsense, nphon*Pinac, nphon*Pcout, Duty, Vimo, Ivao, IIvao.real" self.lvThalfDat_lblstr=["Tphase:Ivao", "Tphase:IIvao.real", "Tphase:Duty"] self.gppPlotter.open_gpp("lvThalfDat") if not useCurrentLoop: Duty=self.get_Duty(self.Vinac) for icnt in range(-nskip, n): self.QCss+=(10.0e-6*self.thalf*(self.get_Vss()<5.0)) if not useCurrentLoop: # assumption-> load current constant iloadavg=self.parent.CL0.QCOUT.get_VCout()/self.parent.Rload ploadavg=self.parent.CL0.QCOUT.get_VCout()*iloadavg Vimo=self.parent.get_Vimo(self.parent.get_Iimo(self.parent.get_Vvinac(self.Vinac), self.Vvao)) #assuming Vvao fluctuations small compared to Vvinac over thalf Iinac=((self.Vvao>1.0)*Vimo/self.parent.Vreg_currsense)*self.parent.get_Iin(sqrt(2)*self.parent.Vac_rms, 2*self.parent.Ptot/2) # assumption-> vinac power commanded by Vimo instantly, and duty known Icout=(1.0-Duty)*(Iinac/Duty)*(float(not self.parent.SkipHalf)) #current from one inductor into cout averaged over dt Icout_net=(1.0+float(self.parent.EnablePh1))*Icout-iloadavg #n inductors minus load current DVCout=np.cumsum(Icout_net*dt)/self.parent.Cout #variation in VCout vs dt VCout=self.parent.CL0.QCOUT.get_VCout()+DVCout self.parent.CL0.QCOUT.QCout += (1.0+float(self.parent.EnablePh1))*self.thalf*np.mean(Icout) - self.thalf*iloadavg #some redundancy, more than one way to calc this Pinac=self.Vinac*Iinac Pcout=VCout*Icout if useCurrentLoop: self.parent.CL0.StatsCL.__init__() self.parent.CL0.stepPWM(0, 1667, plotDatPks=(icnt==n-1), debug=False) phTest, Duty, Vimo, Iinac, Icout, Pinac, Pcout, VCout = self.parent.CL0.StatsCL.getInterpAll(self.Tau) #Iinac,Icout changed def here to instantaneous iloadavg=np.mean(VCout/self.parent.Rload) ploadavg=np.mean(VCout*iloadavg) Ivao=self.get_Ivao(self.get_Vsense(VCout)) IIvao=self.fftIntegrator_Ivao.fftIntegrate(Ivao, flagSavePrev=True) #True:save integrand to history self.Vvao=max(self.Vvao_bounds[0], min(self.Vvao + IIvao[-1].real, self.Vvao_bounds[1])) if plotThalfDat: #realtime display, never to file nphon=(1.0+float(self.parent.EnablePh1)) Tphase, Vsense = self.thalf*self.Tau/np.pi, self.get_Vsense(VCout) self.lvThalfDat=np.vstack(eval("["+self.lvThalfDat_str+"]")) #self.gppPlotter.gpp_Plot(self.lvThalfDat.transpose(), self.lvThalfDat_str, self.lvThalfDat_lblstr, self.gppPlotter.gpp["lvThalfDat"], lineWidth=2.0) self.gppPlotter.gpp_Plot(self.lvThalfDat.transpose(), self.lvThalfDat_str, self.lvThalfDat_lblstr, "lvThalfDat", lineWidth=2.0) if icnt>=0: deltaVao, IvaoAvg, PinacAvg, PcoutAvg, IcoutAvg = IIvao[-1], np.mean(Ivao), np.mean(Pinac), np.mean(Pcout), np.mean(Icout) lvIterDat.append(eval("["+self.lvIterDat_str+"]")) self.addClockThalf() self.parent.SkipHalf=(self.parent.SkipHalfMod>1) and not ((self.cntThalf%self.parent.SkipHalfMod)==0) if plotIterDat: #self.gppPlotter.gpp_Plot(lvIterDat, self.lvIterDat_str, self.lvIterDat_lblstr, lineWidth=1.0, toFile=self.parent.get_FileName("lvIterDat.pdf")) gppconfig="set grid ytics back ; \n set grid lw 0.3 lt 1 lc rgb 0x80888888 ; \n" self.gppPlotter.gpp_Plot(lvIterDat, self.lvIterDat_str, self.lvIterDat_lblstr, lineWidth=1.0, toFile=self.parent.get_FileName("lvIterDat.pdf"), config=gppconfig) ######################################################################################################### class CurrentLoop: Vcaox_bounds=[0.0, 6.0] nptsIcaox=200 #puts 1st stage root finding error < 0.5% Vramp_min=0.7 #not adjustable def __init__(self, parent): self.parent=parent self.fcycIvao=2*self.parent.Fac self.thalf=(1.0/(2*self.parent.Fac)) self.fcycIcaox=self.parent.Fpwm self.tpwm=(1.0/self.parent.Fpwm) self.npwmMax=round(self.parent.Fpwm/(2*self.parent.Fac)) self.ntonPh=np.pi/self.npwmMax self.Tau = np.linspace(0.0, 1.0, self.nptsIcaox, endpoint=False) self.lcDatPks=[] self.lcDatRT=[] self.gppPlotter=gnuPlotter() self.Vimo=None self.StatsCL=self.statsCL() self.QCOUT = self.QCout_Tau_Obj(self, 0.0) self.IL0 = self.ILx_Tau_Obj(self, self.QCOUT) self.CA0 = self.CAx_CSx_Obj(self, self.IL0) self.fftIntegrator_Icaox = fftIntegrator(self.nptsIcaox, self.fcycIcaox, 20, 20, self.get_ZType2_CAOx, self.get_ZType2_CAOx_DC, self.CA0.get_Vcaox, self.Vcaox_bounds) self.Vramp = self.get_Vramp(self.Tau, self.Vramp_min, self.parent.DVramp) self.Clock=self.parent.Clock self.cntTPWM=0 def addClockTPWM(self, n=1): self.cntTPWM+=n self.Clock+=n*1.0/(self.parent.Fpwm) def get_Vramp(self, tau, vmin, dv): return vmin + tau*dv #tau in [0,1] ######################################################################################################### class statsCL: def __init__(self): self.nPh0=[] self.TauOn=[] self.Vimo=[] self.ILx_Vinac=[] self.ILx_VCout=[] self.Pinac=[] self.Pcout=[] self.VCoutAvg=[] def addStats(self, nph0, tauon, vimo, ilxvinac, ilxvcout, pinac, pcout, vcoutavg): self.nPh0.append(nph0) self.TauOn.append(tauon) self.Vimo.append(vimo) self.ILx_Vinac.append(ilxvinac) self.ILx_VCout.append(ilxvcout) self.Pinac.append(pinac) self.Pcout.append(pcout) self.VCoutAvg.append(vcoutavg) def getPCoutAvg(self): return np.mean( np.array(self.VCoutAvg)*np.array(self.ILx_VCout)*(1.0-np.array(self.TauOn)) ) def getInterpAll(self, taunew): r=[ np.interp(taunew, self.nPh0, arrcurr) for arrcurr in [self.nPh0, self.TauOn, self.Vimo, self.ILx_Vinac, self.ILx_VCout, self.Pinac, self.Pcout, self.VCoutAvg] ] return r[0],r[1],r[2],r[3],r[4],r[5],r[6],r[7] ######################################################################################################### class CAx_CSx_Obj: def __init__(self, parentcl, ilxobj, vcaox=0.0): self.parentCL=parentcl self.ILxobj=ilxobj self.Vcsx=None self.Icaox=None self.IIcaox=None self.PWMx=None self.solnFound=False self.TauCurr=None self.TauRoot=None self.TauRootFin=None self.ErrAvg_PWMx=0.0 self.Vcaox=vcaox self.Vcaox_prev=0.0 self.tStamp=-1.0 def reset(self): self.solnFound, self.TauRoot, self.TauRootFin, self.Vcaox_Next = False, None, None, None def get_Vcsx(self, ilx): return ilx*self.parentCL.parent.Rsequiv def get_Icaox(self, vimo, vcsx): return self.parentCL.parent.Gmc*(vimo-vcsx) def get_Vcaox(self): return self.Vcaox def get_PWMx(self): return self.parentCL.Vramp-(self.get_Vcaox() + self.IIcaox.real) def configure(self, nph, tauroot, flagSavePrev): self.ILxobj.configure(nph, tauroot) self.Vcsx=self.get_Vcsx(self.ILxobj.get_ILxCurrArr(self.parentCL.Tau)) self.Icaox=self.get_Icaox(self.parentCL.Vimo, self.Vcsx) self.IIcaox=self.parentCL.fftIntegrator_Icaox.fftIntegrate(self.Icaox, flagSavePrev) self.PWMx=self.parentCL.Vramp-(self.get_Vcaox()+self.IIcaox.real) self.tStamp=self.parentCL.Clock def peakLimitTrigger(self): return (self.Vcsx > self.parentCL.parent.Vreg_currpeak).any() def has_Root(self, arr): return ( (arr > 0.0).any() and (arr < 0.0).any() ) def has_Root_PWMx(self): return self.has_Root(self.PWMx) def has_Root_PeakLimit(self): return self.has_Root(self.Vcsx[:round(self.parentCL.nptsIcaox*self.TauRoot)] - self.parentCL.parent.Vreg_currpeak) def get_FnRoot_PWMxArr(self): return self.PWMx[round(self.parentCL.nptsIcaox*self.TauRoot)] def get_TauRoot_PeakLimit(self, tlo=-0.1, thi=1.1): result=None try: result=findRoot( lambda t: np.interp(t, self.parentCL.Tau, self.Vcsx - self.parentCL.parent.Vreg_currpeak), tlo, thi ) except: pass return result def get_TauRoot_PWMx(self, tlo=-0.1, thi=1.1): result=None try: result = findRoot( lambda t: np.interp(t, self.parentCL.Tau, self.PWMx), tlo, thi ) except: pass return result def get_TauRoot_Any(self, tlo=-0.1, thi=1.1): rootpwm, rootpklim = None, None if self.has_Root_PWMx(): rootpwm = self.get_TauRoot_PWMx(tlo, thi) if self.has_Root_PeakLimit(): rootpklim = self.get_TauRoot_PeakLimit(tlo, thi) elif (rootpwm is None) and self.peakLimitTrigger(): rootpklim = 0.0 if rootpklim is None: #easy, no peaklimit, simple root result = rootpwm, False elif rootpwm is not None: #both roots found, use earliest result = np.fmin(rootpwm, rootpklim), (rootpklim<=rootpwm) #print("Both Roots Found. rootpwm, rootpklim, result == root, is_peak_limited", rootpwm, rootpklim, result) else: #no roots found, returning None or zero result = rootpklim, True print("No Roots Found. rootpwm, rootpklim, result == root, is_peak_limited", rootpwm, rootpklim, result) return result def get_FnRoot_PWMx(self, tauroot): return np.interp(tauroot, self.parentCL.Tau, self.PWMx) def get_IIcaox_Last(self): return self.IIcaox[-1].real def save_Vcaox_Next(self): self.Vcaox_Next=max(self.parentCL.Vcaox_bounds[0], min(self.Vcaox + self.get_IIcaox_Last(), self.parentCL.Vcaox_bounds[1])) def apply_Vcaox_Next(self): self.Vcaox_prev=self.get_Vcaox() self.Vcaox=self.Vcaox_Next ######################################################################################################### class ILx_Tau_Obj: def __init__(self, parentcl, qcobj, ilx=0.0): self.parentCL=parentcl self.Lx=self.parentCL.parent.Lx #convenience copy self.QCobj=qcobj self.ILx=ilx self.Aac=sqrt(2)*self.parentCL.parent.Vac_rms self.Wac=2*np.pi*self.parentCL.parent.Fac self.Kac=self.Aac/(self.Wac*self.Lx) self.ILxCurrArr=None def configure(self, nph0, tauOn): self.nPh0 = nph0 self.TauOn=tauOn self.QCobj.configure(self.TauOn) self.ILxCurrArr=None #must regenerate after config self.set_ILxCurrArr(self.parentCL.Tau) self.ILx_Next=None def get_DILx_Vinac(self, tauI, tauF, vcout=0): #using positive current convention # i(t) ~= -I0 - ( Vinac - VCout )*t/L # Vinac ~= A*sin(ph0 + dph/2); dph=tpwm*TauRoot --> mid interval voltage approx, matters little at small dt # VCout ~= QCout/Cout dphI, dphF = self.parentCL.ntonPh*tauI, self.parentCL.ntonPh*tauF dt=self.parentCL.tpwm*(tauF-tauI) vinac=self.Aac*np.sin(self.nPh0+dphI+dphF/2.0) return (vinac-vcout)*dt/self.Lx def zget_DILx_Vinac(self, tauI, tauF, vcout=0): #using positive current convention #return ( (sqrt(2)*self.parentCL.parent.Vac_rms/(2*np.pi*self.parentCL.parent.Fac))*(np.cos(theta)-np.cos(theta+self.parentCL.ntonPh*tau)) - vcout*self.parentCL.tpwm*tau )/self.Lx # i(t) = -C0 + K*cos(ph0+w*t) + VCout*t/L assuming VCout=0 when t<0; since cos decreasing, i getting more negative as desired # -C0 = -I0 - K*cos(ph0) sets initial current to I0 < 0 dphI=self.parentCL.ntonPh*tauI dphF=self.parentCL.ntonPh*tauF dt=self.parentCL.tpwm*(tauF-tauI) dCos=-(np.cos(self.nPh0+dphF) - np.cos(self.nPh0+dphI)) return self.Kac*dCos - vcout*dt/self.Lx def _get_ILx_On(self, tauF, flagICout=False): if not flagICout: return self.ILx + self.get_DILx_Vinac(0.0, tauF) else: return 0.0 def _get_ILx_Off(self, tauF): result = self._get_ILx_On(self.TauOn) + self.get_DILx_Vinac(self.TauOn, tauF, self.parentCL.QCOUT.get_VCout_Tau(self.TauOn)) return np.fmax(0.0, result) def get_ILx_Tau(self, tau, flagICout=False): #prevent negative currents for root finding, but retain ability to explicitly check for DCM operation result = np.piecewise(tau, \ [tau<0.0, (tau>=0.0)&(tau<=self.TauOn), (tau>self.TauOn)&(tau<=1.0), (tau>1.0)],\ [0, lambda tau: self._get_ILx_On(tau, flagICout), lambda tau: self._get_ILx_Off(tau), 0] ) return np.fmax(result, 0.0) def get_ILxAvgOn(self): return (self.get_ILx_Tau(0.0)+self.get_ILx_Tau(self.TauOn))/2.0 def get_ILxAvgOff(self): return (self.get_ILx_Tau(1.0)+self.get_ILx_Tau(self.TauOn))/2.0 def get_ILx_TauOn(self): return self.get_ILx_Tau(self.TauOn) def get_DE_Vinac(self, tauI, tauF, vcout=0.0): # E(t)-E(0) ~= Vinac*( I0*t + (1/2)*(Vinac - VCout)*t**2/L ) dphI=self.parentCL.ntonPh*tauI dphF=self.parentCL.ntonPh*tauF dt=self.parentCL.tpwm*(tauF-tauI) vinac=self.Aac*np.sin(self.nPh0+dphI+dphF/2.0) i0 = self.get_ILx_Tau(tauI) return vinac*dt*(i0 + 0.5*(vinac-vcout)*dt/self.Lx) def zget_DE_Vinac(self, tauI, tauF, vcout=0.0): # E(t)-E(0) = (A/w)*C(t)*I0 + (1/2)*K*(A/w)*(C(t))**2 - VCout*K*t*C(t) - VCout*(K/w)*S(t) # Don't forget there are 2 inductors! # C(t) = -(cos(ph0+w*t)-cos(ph0)) > 0 # S(t) = (sin(ph0+w*t)-sin(ph0)) # K = (A/(w*L)) dphI=self.parentCL.ntonPh*tauI dphF=self.parentCL.ntonPh*tauF dt=self.parentCL.tpwm*(tauF-tauI) dCos=-(np.cos(self.nPh0+dphF) - np.cos(self.nPh0+dphI)) dSin=(np.sin(self.nPh0+dphF) - np.sin(self.nPh0+dphI)) i0 = self.get_ILx_Tau(tauI) return (self.Aac/self.Wac)*dCos*(i0 + 0.5*self.Kac*dCos) - vcout*self.Kac*(dt*dCos + dSin/self.Wac) def get_DE_VCout(self, tauI, tauF, vcout): # E(t)-E(0) ~= -VCout*( I0*t + (1/2)*(Vinac - VCout)*t**2/L ) dphI=self.parentCL.ntonPh*tauI dphF=self.parentCL.ntonPh*tauF dt=self.parentCL.tpwm*(tauF-tauI) vinac=self.Aac*np.sin(self.nPh0+dphI+dphF/2.0) i0 = self.get_ILx_Tau(tauI) #return -vcout*dt*(i0 + 0.5*(vinac-vcout)*dt/self.Lx) return np.fmin(-vcout*dt*(i0 + 0.5*(vinac-vcout)*dt/self.Lx), 0.0) #DIODE --> ALWAYS NEGATIVE def zget_DE_VCout(self, tauI, tauF, vcout): # E(t)-E(0) = -VCout*(I0 + K*cos(ph0))*t + VCout*(K/w)*S(t) + (1/2)*VCout**2*(t**2)/L dphI=self.parentCL.ntonPh*tauI dphF=self.parentCL.ntonPh*tauF dt=self.parentCL.tpwm*(tauF-tauI) dSin=(np.sin(self.nPh0+dphF) - np.sin(self.nPh0+dphI)) i0 = self.get_ILx_Tau(tauI) return vcout*( (-i0 - self.Kac*np.cos(self.nPh0+dphI) + 0.5*vcout*dt/self.Lx)*dt + self.Kac*dSin/self.Wac ) def set_ILx(self, ilx): self.ILx=max(0.0, ilx) def save_ILx_Next(self): self.ILx_Next=self.get_ILx_Tau(1.0) def apply_ILx_Next(self): tOn=self.parentCL.tpwm*(self.TauOn) deVinacOn=self.get_DE_Vinac(0.0, self.TauOn) tOff=self.parentCL.tpwm*(1.0-self.TauOn) vcoutRoot=self.QCobj.get_VCout_Tau(self.TauOn) deVinacOff=self.get_DE_Vinac(self.TauOn, 1.0, vcoutRoot) self.parentCL.parent.PM0.add_E_T("DE_Vinac_On", deVinacOn, tOn) self.parentCL.parent.PM0.add_E_T("DE_Vinac_Off", deVinacOff, tOff) self.parentCL.parent.PM0.add_E_T("DE_Vinac_Tot", deVinacOn, tOn) self.parentCL.parent.PM0.add_E_T("DE_Vinac_Tot", deVinacOff, tOff) de_ilx_vcout = self.get_DE_VCout(self.TauOn, 1.0, vcoutRoot) self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", deVinacOn, tOn) self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", deVinacOff, tOff) self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", de_ilx_vcout, self.parentCL.tpwm*(1.0-self.TauOn)*0) self.set_ILx(self.ILx_Next) def set_ILxCurrArr(self, tauarr): self.ILxCurrArr = self.get_ILx_Tau(tauarr) def get_ILxCurrArr(self, tauarr): #self.set_ILxCurrArr(tauarr) return self.ILxCurrArr def get_ILxCout_Arr(self): return self.get_ILx_Tau(self.parentCL.Tau, flagICout=True) ######################################################################################################### class QCout_Tau_Obj: def __init__(self, parentcl, qcout=0.0): self.parentCL=parentcl self.QCout=qcout self.Cout=self.parentCL.parent.Cout #convenience copy def configure(self, tauOn): self.TauOn=tauOn self.QCout_Next=None def _get_ICout_On(self, tau): tRC=(self.parentCL.parent.Rload*self.Cout) return -(self.QCout/tRC)*np.exp(-self.parentCL.tpwm*tau/(tRC)) def _get_ICout_Off(self, tau): return self._get_ICout_On(0*tau) def get_ICout_Tau(self, tau, flagPosOk=False): result = np.piecewise(tau, \ [tau<0.0, (tau>=0.0)&(tau<=self.TauOn), (tau>self.TauOn)&(tau<=1.0), (tau>1.0)],\ [0, lambda tau: self._get_ICout_On(tau), lambda tau: self._get_ICout_Off(tau), 0] ) if flagPosOk: return result else: return np.fmin(result, 0.0) def get_DQCout_Rload_dT(self, q0, t): # Q(t)-Q(0) = Q0*( exp(-t/(Rload*Cout)) - 1 ) return q0*(np.exp(-t/(self.parentCL.parent.Rload*self.Cout))-1.0) def get_DE_Rload_dQ(self, qI, qF): # E(t)-E(0) = (1/2)*( Q(t)**2 - Q(0)**2 )/C return 0.5*(qF**2-qI**2)/self.Cout def get_DE_Rload_dT(self, q0, t): # DECout(t) = ELx - t*V(0)**2/Rload ; assumption: constant current return -t*((q0/self.Cout)**2)/self.parentCL.parent.Rload def _get_QCout_On(self, tau): return self.QCout + self.get_DQCout_Rload_dT(self.QCout, self.parentCL.tpwm*tau) def _get_QCout_Off(self, tau): # DECout(t) = ELx - t*V(0)**2/Rload # QCout_new = sqrt( QCout**2 + 2*DECout*Cout ) de_ilx_vcout = -(1.0+float(self.parentCL.parent.EnablePh1))*self.parentCL.IL0.get_DE_VCout(self.TauOn, tau, self._get_QCout_On(self.TauOn)/self.Cout) de_cout_rload = self.get_DE_Rload_dT(self.get_QCout_TauOn(), self.parentCL.tpwm*(tau-self.TauOn)) de_cout_net = de_ilx_vcout + de_cout_rload #notice this does change voltage a bit qnew2=np.fmax(0.0, (self._get_QCout_On(self.TauOn)**2 + 2.0*de_cout_net*self.Cout)) return np.sqrt(qnew2) def get_QCout_Tau(self, tau): #prevent negative charge, ever useful in cap context? result = np.piecewise(tau, \ [tau<0.0, (tau>=0.0)&(tau<=self.TauOn), (tau>self.TauOn)&(tau<=1.0), (tau>1.0)],\ [0, lambda tau: self._get_QCout_On(tau), lambda tau: self._get_QCout_Off(tau), 0] ) return np.fmax(result, 0.0) def get_VCout_Tau(self, tau): return self.get_QCout_Tau(tau)/self.Cout def get_QCout_TauOn(self): return self.get_QCout_Tau(self.TauOn) def get_VCout_TauOn(self): return self.get_QCout_Tau(self.TauOn)/self.Cout def get_PCout_Lx(self): return -self.parentCL.IL0.get_DE_VCout(self.TauOn, 1.0, self.get_VCout_TauOn())/self.parentCL.tpwm #return -self.parentCL.IL0.get_DE_VCout(self.TauOn, 1.0, self._get_QCout_On(self.TauOn)/self.Cout)/self.parentCL.tpwm def get_VCoutAvg(self): return (self.get_VCout_Tau(0.0)+self.get_VCout_Tau(self.TauOn)+self.get_VCout_Tau(1.0))/3.0 def add_QCout_dE(self, dE): # QCout_new = sqrt( QCout**2 + 2*DECout*Cout ) qnew2=self.QCout**2 + 2*dE*self.Cout qnew2=max(0.0, qnew2) self.set_QCout(self.QCout + sqrt(qnew2)) def set_QCout(self, qcout): self.QCout=max(0.0, qcout) def get_VCout(self): return self.QCout/self.Cout def save_QCout_Next(self): self.QCout_Next=self.get_QCout_Tau(1.0) def apply_QCout_Next(self): de_ilx_vcout = -(1.0+float(self.parentCL.parent.EnablePh1))*self.parentCL.IL0.get_DE_VCout(self.TauOn, 1.0, self._get_QCout_On(self.TauOn)/self.Cout) de_cout_rload = self.get_DE_Rload_dT(self.get_QCout_TauOn(), self.parentCL.tpwm*(1.0-self.TauOn)) de_cout_net = de_ilx_vcout + de_cout_rload self.parentCL.parent.PM0.add_E_T("de_ilx_vcout", de_ilx_vcout, self.parentCL.tpwm) self.parentCL.parent.PM0.add_E_T("de_cout_rload", de_cout_rload, self.parentCL.tpwm*(1.0-self.TauOn)) self.parentCL.parent.PM0.add_E_T("de_cout_net", de_cout_net, self.parentCL.tpwm) dq = self.get_DQCout_Rload_dT(self.QCout, self.parentCL.tpwm*self.TauOn) deon = self.get_DE_Rload_dQ(self.QCout, self.QCout+dq) deoff = self.get_DE_Rload_dT(self.get_QCout_TauOn(), self.parentCL.tpwm*(1.0-self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_QCout_Rload_On", deon, self.parentCL.tpwm*(self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_QCout_Rload_Off", deoff, self.parentCL.tpwm*(1.0-self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_QCout_Rload_Tot", deon+deoff, self.parentCL.tpwm) self.set_QCout(self.QCout_Next) #----------------------------------------------------------------------------------------------------------------------------------- def stepPWM(self, n0=None, n=1, plotDatPks=False, plotDatTau=False, plotDatRT=False, debug=False, dtau=1.0/(2.0*nptsIcaox), eps=1.0e-6): self.gppPlotter.open_gpp("lcDatTau") self.lcDatTau=np.array([]).reshape(7,0) self.ILxCurr1, self.ILxCout1, self.ILxCurrTot, self.ILxCoutTot = self.Tau, self.Tau, self.Tau, self.Tau, self.lcDatTau_str= "self.Clock+self.tpwm*self.Tau, self.IL0.ILxCurrArr, self.IL0.get_ILxCout_Arr(), self.ILxCurr1, self.ILxCout1, self.ILxCurrTot, self.ILxCoutTot" self.lcDatTau_lblstr=["Tau:ILxCurrArr, Tau:ILxCurr1, Tau:ILxCurrTot","Tau:get_ILxCout_Arr(), Tau:ILxCout1, Tau:ILxCoutTot"] self.gppPlotter.open_gpp("lcDatRT") self.lcDatRT_str= "self.Tau, self.IL0.ILxCurrArr, self.CA0.Vcsx, self.Vimo, self.CA0.Icaox, self.CA0.IIcaox.real, self.CA0.PWMx" self.lcDatRT_lblstr=["Tau:ILxCurrArr","Tau:PWMx"] self.gppPlotter.open_gpp("lcDatPks") self.lcDatPks_str= "self.Clock, self.IL0.get_ILx_Tau(0.0), self.IL0.get_ILx_TauOn(), self.IL0.get_ILx_Tau(1.0),"+ \ "self.QCOUT.get_VCout_Tau(0.0), self.QCOUT.get_VCout_TauOn(), self.QCOUT.get_VCout_Tau(1.0),"+ \ "self.CA0.Vcaox_prev, self.CA0.get_Vcaox(), self.CA0.get_IIcaox_Last(), self.CA0.Icaox[-1],"+ \ "ncnt, self.CA0.get_FnRoot_PWMxArr(), self.CA0.get_FnRoot_PWMx(self.CA0.TauRootFin), self.CA0.TauRoot, 2*self.CA0.TauRootFin, self.CA0.ErrAvg_PWMx, "#+ \ self.lcDatPks_lblstr=["Clock:VCout_TauOn(), Clock:VCout_Tau(1.0)", "Clock:ILx_Tau(0.0), Clock:ILx_TauOn(), Clock:ILx_Tau(1.0), Clock:TauRootFin"] if n0 is not None: self.cntTPWM=n0 if n0==0: self.parent.PM0.__init__() for ncnt in range(self.cntTPWM, self.cntTPWM+n): self.CA0.reset() nPh0 = (ncnt%self.npwmMax)*self.ntonPh #this + delta must be <= Pi; recall pain around integral of fabs(sin) at Pi Vinac=self.parent.get_Vinac_Sin(nPh0 + self.ntonPh*self.Tau) self.Vimo=self.parent.get_Vimo(self.parent.get_Iimo(self.parent.get_Vvinac(Vinac), self.parent.VL0.Vvao)) #assuming Vvao fluctuations small compared to Vvinac over tpwm if (self.CA0.get_Vcaox()>self.Vramp_min) and (not self.parent.SkipHalf): #begins with state On if debug: print(" ON, trying Dmax for root find; ") self.CA0.TauRoot=self.parent.Dmax else: #off entire cycle if debug: print(" OFF entire cycle. Accepting 0 solution as final... ILx:", self.IL0.ILx, self.QCOUT.QCout) self.CA0.TauRoot = 0.0 #WOW USING INTEGER 0 HERE BROKE IT IN SNEAKIEST SHITTY WAY, not impressed, python FAHHHHKKK self.CA0.solnFound = True self.CA0.configure(nPh0, self.CA0.TauRoot, flagSavePrev=False) #False:not committed to this integrand yet, don't save in history if not self.CA0.solnFound: if self.CA0.has_Root_PWMx(): self.CA0.TauRoot, self.CA0.PeakLimit = self.CA0.get_TauRoot_Any() if self.CA0.PeakLimit: #print("PEAK LIMIT:", ncnt, self.CA0.TauRoot) self.CA0.solnFound=True self.CA0.TauRootFin=self.CA0.TauRoot elif (self.CA0.TauRoot < self.parent.Dmax) and (self.CA0.TauRoot > 0.0): #Good Root Found, refine and calc actual solution self.CA0.configure(nPh0, self.CA0.TauRoot, flagSavePrev=False) self.CA0.TauRootFin = self.CA0.get_TauRoot_PWMx(self.CA0.TauRoot-5*dtau, self.CA0.TauRoot+5*dtau) if ( fabs(self.CA0.get_FnRoot_PWMx(self.CA0.TauRootFin)) < min(fabs(self.CA0.get_FnRoot_PWMxArr()),eps) ) and ( fabs(fabs(self.CA0.TauRootFin) - fabs(self.CA0.TauRoot)) < dtau ): #accept new tau_root pwmx_root=fabs(self.CA0.get_FnRoot_PWMx(self.CA0.TauRootFin)) else: print("ERROR?? Should Never Happen. keeping old tau_root: pwmx_root_ini={:1.6e}, pwmx_root_fin={:1.6e}, tau_root={:1.6e}, tau_root_fin={:1.6e}".format(self.CA0.get_FnRoot_PWMxArr(), self.CA0.get_FnRoot_PWMx(self.CA0.TauRootFin), self.CA0.TauRoot, self.CA0.TauRootFin)) self.CA0.TauRootFin=None pwmx_root=fabs(self.CA0.get_FnRoot_PWMxArr()) quit() self.CA0.ErrAvg_PWMx=(self.CA0.ErrAvg_PWMx + pwmx_root)/2.0 self.CA0.solnFound = True if not self.CA0.solnFound: if self.CA0.peakLimitTrigger(): self.CA0.TauRoot=0.0 else: self.CA0.TauRoot=self.parent.Dmax if debug or 1: if self.CA0.TauRoot > self.parent.Dmax: print(" OOB: Root exceeds bounds: Dmax limiting; tau_root=", ncnt, self.CA0.TauRoot, end=" \n") else: pass #print(" No Root found in PWMx: peak limited-->zero; tau_root=", ncnt, self.CA0.TauRoot, end=" \n") self.CA0.solnFound = True if self.CA0.solnFound: if self.CA0.TauRootFin is None: self.CA0.TauRootFin=self.CA0.TauRoot self.CA0.configure(nPh0, self.CA0.TauRootFin, flagSavePrev=True) #self.CA0.update_Vcaox_with_IIcaox() self.StatsCL.addStats(nPh0, self.CA0.TauRootFin, np.mean(self.Vimo), self.IL0.get_ILxAvgOn(), self.IL0.get_ILxAvgOff(), np.mean(Vinac)*self.IL0.get_ILxAvgOn(), self.QCOUT.get_PCout_Lx(), self.QCOUT.get_VCoutAvg()) #WHY HERE? some variables lose scope outside loop self.lcDatPks.append(eval("["+self.lcDatPks_str+"]")) if plotDatTau: self.lcDatTau=np.hstack([self.lcDatTau, np.vstack(eval("["+self.lcDatTau_str+"]"))]) else: print("ERROR: no root found or assumed. stopping") quit() #since some are dependent, must save all new values first self.CA0.save_Vcaox_Next() self.IL0.save_ILx_Next() self.QCOUT.save_QCout_Next() #now ok to apply changes self.CA0.apply_Vcaox_Next() self.IL0.apply_ILx_Next() #also updates power meter self.QCOUT.apply_QCout_Next() #also updates power meter self.addClockTPWM() if plotDatRT: #realtime display, never to file #self.lcDatRT=list( eval("zip(*["+self.lcDatRT_str+"])") ) self.lcDatRT=np.vstack(eval("["+self.lcDatRT_str+"]")) #self.gppPlotter.gpp_Plot(self.lcDatRT.transpose(), self.lcDatRT_str, self.lcDatRT_lblstr, self.gppPlotter.gpp["lcDatRT"]) self.gppPlotter.gpp_Plot(self.lcDatRT.transpose(), self.lcDatRT_str, self.lcDatRT_lblstr, "lcDatRT") if plotDatTau: #contortions needed to dup currents pi out of phase and add them self.lcDatTau[3]=np.interp(self.lcDatTau[0]-self.tpwm/2.0, self.lcDatTau[0], self.lcDatTau[1]) self.lcDatTau[4]=np.interp(self.lcDatTau[0]-self.tpwm/2.0, self.lcDatTau[0], self.lcDatTau[2]) self.lcDatTau[5]=self.lcDatTau[1]+self.lcDatTau[3] self.lcDatTau[6]=self.lcDatTau[2]+self.lcDatTau[4] tmpname=self.parent.get_FileName("ilxTau_{}_{}.pdf".format(n, self.cntTPWM)) #self.gppPlotter.gpp_Plot(self.lcDatTau.transpose(), self.lcDatTau_str, self.lcDatTau_lblstr, self.gppPlotter.gpp["lcDatTau"], lineWidth=0.5, toFile=tmpname) self.gppPlotter.gpp_Plot(self.lcDatTau.transpose(), self.lcDatTau_str, self.lcDatTau_lblstr, "lcDatTau", lineWidth=0.5, toFile=tmpname) if plotDatPks: #self.gppPlotter.gpp_Plot(self.lcDatPks, self.lcDatPks_str, self.lcDatPks_lblstr, self.gppPlotter.gpp["lcDatPks"], lineWidth=1.0, toFile=self.parent.get_FileName("lcDatPks.pdf")) self.gppPlotter.gpp_Plot(self.lcDatPks, self.lcDatPks_str, self.lcDatPks_lblstr, "lcDatPks", lineWidth=1.0, toFile=self.parent.get_FileName("lcDatPks.pdf")) self.parent.PM0.Report() return self.lcDatPks[-1] def get_ZType2_CAOx(self, bk, f): return self.parent.get_Zcomp_type2(bk, f, self.parent.Rzc, self.parent.Czc, self.parent.Cpc) def get_ZType2_CAOx_DC(self): # should be integration coefficient, simply 1/(cz+cp), to calc DC component integral approx over dt: dv = dq/Ctot ~= dt*DCcomponentOfCurrent/(cp+cz) return self.parent.get_Zcomp_type2(1, 1/(2*pi*1j), 0, self.parent.Czc, self.parent.Cpc) ######################################################################################################### # FILTERS : Complex Impedance, Bode Plots def get_Zcomp_type2(self, bk, f, rz, cz, cp): omega = 2*pi*f zz = rz + 1.0/(omega*cz*1j) return bk / ( (omega*cp*1j) + (1.0/zz) ) def get_Z_lst(self, f, res): return [f, res, abs(res), atan(res.imag/res.real)] def get_impedance_div2R2C(self, omega, ra, rb, ca, cb): tmpa = 1.0 / ( (omega*ca*1j) + (1.0/ra) ) tmpb = 1.0 / ( (omega*cb*1j) + (1.0/rb) ) res = tmpb / (tmpa + tmpb) return [omega, res, abs(res), atan2(res.imag, res.real)] ################################## # main ########################### ################################## def main(): # if(len(sys.argv)>1): # l=sys.argv.copy() # l.remove(l[0]) pfc0=UCC28070(plotpdf=True) def formatCustom(a): if type(a) is float: return "{:20.6}".format(a) elif (type(a) is int) or (type(a) is bool): return "{:20}".format(a) else: return "{}".format(a) with open("logs/log_Vac{}_Ptot{}_Lx{}_Fpwm{}_RBac{}_Cout{}.dat".format(\ pfc0.Vac_rms, pfc0.Ptot, round(pfc0.Lx*1.0e6), pfc0.Fpwm, round(pfc0.RBac), round(pfc0.Cout*1.0e6) ), 'w') as logfile: for k in pfc0.__dict__.keys(): #print(k,":",pfc0.__dict__[k], file=logfile) print("{:20} : {}".format(k, formatCustom(pfc0.__dict__[k])), file=logfile) print("ICout at Pmax:", pfc0.get_ICout(sqrt(2)*pfc0.Vac_rms, 2*pfc0.Ptot), "Vripple:", pfc0.get_Vripple(), file=logfile) print("Vvinac max:", pfc0.get_Vvinac(sqrt(2)*pfc0.Vac_rms), "kvff:", pfc0._get_Kvff(pfc0.get_Vvinac(sqrt(2)*pfc0.Vac_rms)), file=logfile) print("get_Vinac_Correction_Factor():", pfc0.get_Vinac_Correction_Factor(), file=logfile) print("IMO_max:", pfc0.get_IMO_max(), "Vimo_max:", pfc0.Rimo*pfc0.get_IMO_max(), file=logfile) print("Duty at top:", pfc0.get_Duty(sqrt(2)*pfc0.Vac_rms), " at peak voltage:", sqrt(2)*pfc0.Vac_rms, file=logfile) print("ILxTop:", pfc0.get_ILxTop(pfc0.Vac_rms,pfc0.Ptot/2), "dILxTop:", pfc0.get_dILxTop(pfc0.Vac_rms), "ILxPeak:", pfc0.get_ILxPeak(), file=logfile) print("dILxTop/ILxTop:", pfc0.get_dILxTop(pfc0.Vac_rms)/pfc0.get_ILxTop(pfc0.Vac_rms,pfc0.Ptot/2), file=logfile) print("Rsequiv:", pfc0.Rsequiv, " --> need Vreg_currpeak >= ", pfc0.Vreg_currsense*pfc0.get_ILxPeak()/pfc0.get_ILxTop(pfc0.Vac_rms, pfc0.Ptot/2), file=logfile) print("LxMin_dILx:", pfc0.get_LxMin_dILx(pfc0.Vac_rms, pfc0.Ptot/2, pfc0.Fpwm, 0.3), "for fracilx=0.3", file=logfile) print("LxMin_CCM at Pmax/4:", pfc0.get_LxMin_CCM(pfc0.Vac_rms, pfc0.Ptot/(2*4), pfc0.Fpwm), file=logfile) print("LxMin_CCM at Pmax/3:", pfc0.get_LxMin_CCM(pfc0.Vac_rms, pfc0.Ptot/(2*3), pfc0.Fpwm), file=logfile) print("PoutMin_CCM:", pfc0.get_PoutMin_CCM(pfc0.Vac_rms, pfc0.Fpwm), "with Lx:", pfc0.Lx, file=logfile) #print("_get_ILx_Vals at n=833:", pfc0._get_ILx_Vals(833, pfc0.Ptot/2.0), file=logfile) #quit() pdiv=1 pfc0.Rload=225.0*pdiv pfc0.EnablePh1=True pfc0.SkipHalfMod=1 pfc0.Css=pfc0.get_Css(0.3) pfc0.CL0.QCOUT.QCout=340*pfc0.CL0.QCOUT.Cout #pfc0.CL0.QCOUT.QCout=600*pfc0.CL0.QCOUT.Cout pfc0.VL0.QCss=1.8*pfc0.Css #pfc0.VL0.QCss=4.7*pfc0.Css pfc0.VL0.Vvao=1.0 + (float(pfc0.SkipHalfMod))*pfc0.DVvao/pdiv pfc0.CL0.CA0.Vcaox=4.6/3 pfc0.CL0.IL0.ILx=0.0 iIn, iFin = 0, round(1667/2) if 0: iStep=33 for ii in [i for i in range(iIn, int(iFin), iStep)]: print("ii:",ii) pfc0.CL0.stepPWM(n=iStep, plotDatPks=0, plotDatTau=1, debug=False) quit() if 1: pfc0.VL0.stepTHALF(60, useCurrentLoop=True) quit() if 0: pfc0.CL0.stepPWM(iIn, iFin-iIn, plotDatPks=1, plotDatTau=1, plotDatRT=0, debug=False) quit() pfc0.plotBode_VoltageLoop() pfc0.plotBode_CurrentLoop() ################################## ################################## ################################## if __name__ == '__main__': main() ###################################################### # Voltage Compensation Network Gain: Gcv(f): # #series Rz Cz: # Zser = Rz + 1/(j*2*pi*f*Cz) # Zser = ( Rz*j*2*pi*f*Cz + 1 )/(j*2*pi*f*Cz) #parallel Zs Cp: # Zpar = 1 / ( 1/Zser + (j*2*pi*f*Cp) ) # Zpar = 1 / ( (j*2*pi*f*Cz) / ( Rz*j*2*pi*f*Cz + 1 ) + (j*2*pi*f*Cp) ) # Zpar = ( Rz*j*2*pi*f*Cz + 1 ) / ( (j*2*pi*f*Cz) + (j*2*pi*f*Cp) * ( Rz*j*2*pi*f*Cz + 1 ) ) # Zpar = ( Rz*j*2*pi*f*Cz + 1 ) / ( (j*2*pi*f*Cz) + (j*2*pi*f*Cp) + (j*2*pi*f*Cp)*(Rz*j*2*pi*f*Cz) ) # Zpar = ( Rz*j*2*pi*f*Cz + 1 ) / ( (j*2*pi*f)*(Cz+Cp) + (j*2*pi*f)*(Cz+Cp)*(j*2*pi*f*Cz*Cp*Rz)/(Cz+Cp) ) # Zpar = ( Rz*j*2*pi*f*Cz + 1 ) / ( (j*2*pi*f)*(Cz+Cp) * ( 1 + (j*2*pi*f*Cz*Cp*Rz)/(Cz+Cp) ) # # Gcv(f) = dVvao/dVout = Kdivvs * Gmv * Zpar # ###################################################### # Voltage Loop Power Stage Gain: Gpsv(f): # # Gpsv(f) = (Ptot/Vout)*(1/(j*2*pi*f*Cout))/DVvao # ###################################################### # # # # # ######################################################################################################## ######################################################################################################## # ENERGY BASED APPROACH TO FIND FAST PWM TIMESCALE INDUCTOR CURRENT VALUES # control strategy: average current must follow # Iin(t)=A*sin(omega*t) where A=sqrt(2)*(Prms/Vrms), omega*t in [0,pi] # the average current must track the sin wave, # ######################################################################################################## ######################################################################################################## # when general, v(t) is A*sin(w*t) - VCout, so # integral di/dt = integral -(A*sin(w*t)-VCout)/L # --> i(t) = -C + ((A/w)*cos(w*t)+VCout*t)/L # but to satisfy initial condition I(0) = -I0, # i(0) = -I0 = -C + ((A/w)*cos(ph0)+VCout*(0)/L) ; assuming VCout off until t=0; # -I0 = -C + ((A/w)*cos(ph0))/L # -C = -I0 - ((A/w)*cos(ph0))/L # at time t later: # i(t) = -C + ((A/w)*cos(ph0+w*t)+VCout*t)/L # i(t) = -I0 - ((A/w)*cos(ph0))/L + ((A/w)*cos(ph0+w*t) + VCout*t)/L ######################################################################################################## general form, exact for all t: # i(t) = -I0 - (A/(w*L))*(cos(ph0) - cos(ph0+w*t)) + VCout*t/L ####################################################################################################################################### # SIMPLEST: for t<<1/Wac # i(t) ~= -I0 - ( Vinac - VCout )*t/L # Vinac ~= A*sin(ph0 + dph/2); dph=tpwm*TauRoot --> mid interval voltage approx # VCout ~= QCout/Cout ####################################################################################################################################### # since cos is decreasing, the current is more negative (into the terminal) from positive ac voltage applied, OK. # # interpreting signs: the current from the positive AC terminal is INTO the terminal, hence negative; # this current flows through the inductor, out of the VCout terminal, doing POSITIVE work on Cout. OK. # the positive VCout on the opposite side changes the current in the opposite direction, OK. # at left terminal, current flowing in against the positive potential does positive work on the inductor, # at right terminal, current flowing out with the positive potential does negative work on the inductor (doing positive work on capacitor). # #------------------------------------------------------------------------------------------------------------------------------ # NOTE CONVENTION: in software, sign convention is opposite of above, ILx is POSITIVE always. # +VCout with positive current inductor-->cap does positive work # consistent with positive AC voltage applied to opposite side of inductor, increasing current in positive direction #------------------------------------------------------------------------------------------------------------------------------ # # Integrate I*V to get energy - simple right? just watch # expand current to convenient form, # i(t) = -C0 + K*cos(ph0+w*t) + VCout*t/L ; # K = (A/(w*L)) # -C0 = -I0 - K*cos(ph0) # consider work integrals on left,right separately: # v_left(t) = A*sin(ph0+w*t) v_right(t)= VCout # LEFT SIDE: positive work done on inductor --> E(t) = - integral v_left(t) * i(t) # E(t) = -integral ( A*sin(ph0+w*t) ) * ( -C0 + K*cos(ph0+w*t) + VCout*t/L )dt # E(t) = -integral ( -A*sin(ph0+w*t)*C0 + A*sin(ph0+w*t)*K*cos(ph0+w*t) + A*sin(ph0+w*t)*VCout*t/L )dt # E(t) = - ( (A/w)*cos(ph0+w*t)*C0 - (1/2)*(A/w)*cos**2(ph0+w*t)*K + ( integral t*sin(ph0+w*t)*A*VCout/L dt ) ) ; #TRICKY, cos**2 not -sin**2; integration constants gonna getchya # by parts: integral t*sin(ph0+w*t) dt = -(1/w)*t*cos(ph0+w*t) - integral -(1/w)*cos(ph0+w*t)*1 dt # by parts: integral t*sin(ph0+w*t) dt = -(1/w)*t*cos(ph0+w*t) + (1/w**2)*sin(ph0+w*t) # E(t) = - ( (A/w)*cos(ph0+w*t)*C0 - (1/2)*(A/w)*cos**2(ph0+w*t)*K + (A*VCout/L)*( -(1/w)*t*cos(ph0+w*t) + (1/w**2)*sin(ph0+w*t) ) ) # E(t) = -(A/w)*cos(ph0+w*t)*C0 + (1/2)*(A/w)*cos**2(ph0+w*t)*K + VCout*(A/(w*L))*t*cos(ph0+w*t) - VCout*(1/w)*(A/(w*L))*sin(ph0+w*t) # E(t)-E(0) = -(A/w)*cos(ph0+w*t)*C0 + (1/2)*(A/w)*cos**2(ph0+w*t)*K + VCout*K*t*cos(ph0+w*t) - VCout*(1/w)*K*sin(ph0+w*t) #MUST EVAL E(t)-E(0); # E(t)-E(0) = (A/w)*C(t)*C0 + (1/2)*(A/w)*(cos**2(ph0+w*t)-cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t) # C(t) = -(cos(ph0+w*t)-cos(ph0)) > 0, S(t) = (sin(ph0+w*t)-sin(ph0)) # E(t)-E(0) = (A/w)*C(t)*(I0+K*cos(ph0)) + (1/2)*(A/w)*(cos**2(ph0+w*t)-cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t) # -C0 = -I0 - K*cos(ph0) # E(t)-E(0) = (A/w)*C(t)*I0 - (A/w)*(cos(ph0+w*t)-cos(ph0))*K*cos(ph0) + (1/2)*(A/w)*(cos**2(ph0+w*t)-cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t) # E(t)-E(0) = (A/w)*C(t)*I0 - (A/w)*(cos(ph0+w*t)*cos(ph0) - cos**2(ph0))*K + (1/2)*(A/w)*(cos**2(ph0+w*t) - cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t) # E(t)-E(0) = (A/w)*C(t)*I0 - (A/w)*(cos(ph0+w*t)*cos(ph0))*K + (1/2)*(A/w)*(cos**2(ph0+w*t) + cos**2(ph0))*K - VCout*K*t*C(t) - VCout*(1/w)*K*S(t) # E(t)-E(0) = (A/w)*C(t)*I0 + (1/2)*K*(A/w)*( -2*cos(ph0+w*t)*cos(ph0) + cos**2(ph0+w*t) + cos**2(ph0) ) - VCout*K*t*C(t) - VCout*(1/w)*K*S(t) # E(t)-E(0) = (A/w)*C(t)*I0 + (1/2)*K*(A/w)*(cos(ph0) - cos(ph0+w*t))**2 - VCout*K*t*C(t) - VCout*(1/w)*K*S(t) ######################################################################################################## # WORK ON INDUCTOR FROM VINAC (these need review - not sure if individual term signs are correct?) ######################################################################################################## # E(t)-E(0) = (A/w)*C(t)*I0 + (1/2)*K*(A/w)*(C(t))**2 - VCout*K*t*C(t) - VCout*(K/w)*S(t) # Don't forget there are 2 inductors! # C(t) = -(cos(ph0+w*t)-cos(ph0)) > 0 # S(t) = (sin(ph0+w*t)-sin(ph0)) # K = (A/(w*L)) # -C0 = -I0 - K*cos(ph0) sets initial current to I0 < 0 # i(t) = -C0 + K*cos(ph0+w*t) + VCout*t/L assuming VCout=0 when t<0; since cos decreasing, i getting more negative as desired ####################################################################################################################################### # SIMPLEST: for t<<1/Wac # E(t)-E(0) ~= Vinac*( I0*t + (1/2)*(Vinac - VCout)*t**2/L ) # Vinac ~= A*sin(ph0 + dph/2); dph=tpwm*TauRoot --> mid interval voltage approx # VCout ~= QCout/Cout ####################################################################################################################################### # E(t)-E(0) looks sensible, positive definite when VCout=0, tends to 0 as t-->0, OK. # continue, # RIGHT SIDE: negative work done by (by definition negative) current flowing out of inductor into +VCout --> E(t) = integral v_right(t) * i(t) # with v_right(t) = VCout when transistor is off (zero otherwise) # E(t) = integral ( VCout ) * ( -C0 + K*cos(ph0+w*t) + VCout*t/L )dt # E(t) = VCout * integral ( -C0 + K*cos(ph0+w*t) + VCout*t/L )dt # E(t) = VCout * ( -C0*t + (1/w)*K*sin(ph0+w*t) + (1/2)*VCout*t**2/L ) #MUST EVAL E(t)-E(0); # E(t)-E(0) = VCout * ( -C0*t + (1/w)*K*S(t) + (1/2)*VCout*(t**2)/L ) # E(t)-E(0) = -VCout*C0*t + VCout*(K/w)*S(t) + (1/2)*VCout**2*(t**2)/L ######################################################################################################## # WORK ON INDUCTOR FROM VCOUT ######################################################################################################## # E(t)-E(0) = -VCout*(I0 + K*cos(ph0))*t + VCout*(K/w)*S(t) + (1/2)*VCout**2*(t**2)/L ####################################################################################################################################### # SIMPLEST: for t<<1/Wac # E(t)-E(0) ~= -VCout*( I0*t + (1/2)*(Vinac - VCout)*t**2/L ) # Vinac ~= A*sin(ph0 + dph/2); dph=tpwm*TauRoot --> mid interval voltage approx # VCout ~= QCout/Cout ####################################################################################################################################### # Note this minus the energy consumed by the load is the energy left to charge the capacitor. # negative in initial current, good. # linear decrease in current yields +quadratic term. # work done by Vinac also adds a +S(t) term; # note this term cancels when two work expressions are added to get total work on inductor # # WORK ON CAPACITOR FROM INDUCTOR # is simply the negative of the previous expression MINUS the load energy demand during that period # # WORK ON LOAD FROM VCOUT # while the inductor is grounded with Vinac working on it, # the load is simply connected to VCout and during this time can be simply modelled as # an exponential decay solution with time constant Rload*Cout, since for Cout # dQ/dt = - V(t)/Rload ; V(t) = Q(t)/Cout # dQ/dt = - Q(t)/(Rload*Cout) --> # Q(t) = Q0*exp(-t/(Rload*Cout)) ; Q0 = Q(0) = V(0)*Cout ######################################################################################################## # Q(t)-Q(0) = Q0*( exp(-t/(Rload*Cout)) - 1 ) ######################################################################################################## # SIMPLEST: no need to change the following, already good. ######################################################################################################## # and the work done by the cap on the load during this time changes the cap energy as # E(t)-E(0) = (1/2)*( Q(t)**2 - Q(0)**2 )/C # # While the inductor is connected to the cap, quickly dumping energy into it, # the load continues to demand energy at a rate # Pload = V(t)**2/Rload # diverting some of the energy of the inductor away from the capacitor. ##### Since the inductor current switching time is about 1e5 times faster than Cout/Rload timescale, ##### it may be ok to assume constant voltage and do energy jump calculations, ##### rather than develop full complicated power flow diffyQ calculations. so, # for a known time t, we can calculate the energy demand as # Eload = Pload * t = t*V(t)**2/Rload # and knowing how much energy ELx the inductor is transferring from above, we find # the energy left for the cap as # DECout = ELx - Eload # DECout(t) = ELx - t*V(t)**2/Rload # which may be positive or negative. If V(t) is changing only slightly, we can use V(0) = Q(0)/Cout, ######################################################## Don't forget 2 inductors contribute! # DECout(t) = ELx - t*V(0)**2/Rload ######################################################## # to update QCout via ECout_new = ECout + DECout, # (1/2)*QCout_new**2/Cout = (1/2)*QCout**2/Cout + DECout # QCout_new**2 = QCout**2 + 2*DECout*Cout ######################################################## # QCout_new = sqrt( QCout**2 + 2*DECout*Cout ) ######################################################## # we can use the new Q to revise the V(t) we used, or just let it be an artificial inefficiency. # we will be slightly underestimating the voltage and the energy increase during inductor dump. # approx 1e-5 fractional error in charge, 10s of millivolts or less. # ######################################################################################################## ######################################################################################################## # transconductance weirdness # like a resistor in series, with current defined by voltage: # iout = C * vin ; C units 1/r # r(vin) = vin/iout # but if iout into non zero va, # r(vin) = (vin-va)/iout = (vin-va)/(C*vin) # r(vin) = (1-(va/vin))/C # which becomes like negative resistance for va>vin #-----------------------------------------------------------------------------------------------------------------------------------------------