#!/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 ######################################################################################################### class gnuPlotter: def __init__(self): self.gpp={} def open_gpp(self, gppname=None): if (gppname is None) or (gppname 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 gppname is None: return gppnew else: self.gpp[gppname]=gppnew def gpp_Plot(self, ldat, strdat, lstrlbl, gppobj=None): 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))] return plotStr+", ".join(litemStr)+" \n" lplotStr=[get_plotStr(lbl) for lbl in lstrlbl] if gppobj is None: gppCurr=self.open_gpp() else: gppCurr=gppobj gppCurr.stdin.write("$idat << EOD \n") for irow in ldat: gppCurr.stdin.write(("{} "*len(irow)+" \n").format(*irow)) gppCurr.stdin.write("EOD \n") 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]+" pause -1 \n") gppCurr.stdin.write("unset multiplot \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): ####################### # 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; start with Vripple instead? self.tVdrop=(1/2)/self.Fac # Power Total W RMS self.Ptot=1600 self.Eff=1.0 # Frequency Switching PWM self.Fpwm=60*3334 #use int multiple of Fac self.Fsync=2*self.Fpwm self.Dmax=0.98 self.fracILx=0.3 #allowed variation of ILx as fraction of ILxAverage; ILxRange=ILx*(1 \pm fracILx/2) self.Enable_Synthesis=0 #disable synthesis if full-time current sensing used ####################### # INTERNAL PARAMETERS ####################### # Regulation Voltages self.Vreg_vsense=3.0 #not adjustable self.DVvao=3.8 #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.3 #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.8e-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 scale vinac to full 3V range, most useful for 120VAC applications!! more accurately detect zero crossing and better signal/noise!! hooray for no current synthesis) self.RAac=3.0e6; self.RBac=self.get_RB(self.RAac, self.Vreg_vsense, self.Vout) 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() # 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.Rsequiv=self.get_Rsequiv() 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.Rrt, self.Rdmx = self.get_Rrt_Rdmx() # Vout filter self.Cout=self.get_Cout(self.Vdrop, self.tVdrop) # VAO compensation type 2 self.fracDVvao=0.02 #fraction of DVvao to scale 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 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.Vvao=0 self.QCss=0 self.Clock=0 self.cntThalf=0 self.Rload=225.0 self.LoopDatV=[] self.nptsIvao=20 self.fcycIvao=2*self.Fac self.gpp={} 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.CL0 = self.CurrentLoop(self) self.PM0 = self.PowerMeter() ######################################################################################################### # OSCILLATOR : Rate, Duty Max def get_Rrt_Rdmx(self, tsyncext=0): rrt = (15.0e9 / self.Fsync) rdmx = rrt * (2.0*self.Dmax - 1.0 - (tsyncext*self.Fsync)) 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_Kdivvs(self): return self.RBvs/(self.RAvs+self.RBvs) def get_Vsense(self, vout): return vout*self.get_Kdivvs() 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.Gmv-iscp, lambda dum: dum*self.Gmv, isat] \ ) return rnew 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 def get_Vss(self): return self.QCss/self.Css def get_Vvao(self): return self.Vvao def get_Cout(self, vdrop, dt): #valid for 'small' vdrop << Vout return (dt*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_Kdivvs()*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_Kdivvs()*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_Kdivvs() * 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) #----------------------------------------------------------------------------------------------------------------------------------- def iter_Loop_Voltage(self, n, useCurrentLoop, nskip=0): if(self.Vvao<=1.0): self.Vvao=1.0001 iterdat=[] if not "iterLV" in self.gpp.keys(): self.gpp["iterLV"]=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE, stderr=subprocess.PIPE) self.gpp["iterLV"].stdin.write("set term x11 size 1900,950; \n") thalf=(1.0/self.fcycIvao) dt=thalf*(1.0/self.nptsIvao) Tau = np.linspace(0, np.pi, self.nptsIvao, endpoint=False) Vinac=self.get_Vinac_Sin(Tau) if not useCurrentLoop: Duty=self.get_Duty(Vinac) for icnt in range(-nskip, n): self.QCss+=(10.0e-6*thalf*(self.get_Vss()<5.0)) #self.Rload=225*(1.5-0.5*sin(2*pi*icnt/14)) if not useCurrentLoop: # assumption-> load current constant iloadavg=self.CL0.QCOUT.get_VCout()/self.Rload ploadavg=self.CL0.QCOUT.get_VCout()*iloadavg Vimo=self.get_Vimo(self.get_Iimo(self.get_Vvinac(Vinac), self.Vvao)) #assuming Vvao fluctuations small compared to Vvinac over thalf # assumption-> vinac power commanded by Vimo instantly, and duty known Iinac=((self.Vvao>1.0)*Vimo/self.Vreg_currsense)*self.get_Iin(sqrt(2)*self.Vac_rms, 2*self.Ptot/2) Icout=(1.0-Duty)*(Iinac/Duty) #current from one inductor into cout averaged over dt Icout_net=2*Icout-iloadavg #2 inductors minus load current DVCout=np.cumsum(Icout_net*dt)/self.Cout #variation in VCout vs dt VCout=self.CL0.QCOUT.get_VCout()+DVCout self.CL0.QCOUT.QCout += 2*thalf*np.mean(Icout) - thalf*iloadavg #some redundancy, more than one way to calc this Pinac=Vinac*Iinac Pcout=VCout*Icout if useCurrentLoop: self.CL0.StatsCL.__init__() self.CL0.stepPWM(0, 1667, plotdat=(icnt==n-1), debug=False) phTest, Duty, Vimo, Iinac, Icout, Pinac, Pcout, VCout = self.CL0.StatsCL.getInterpAll(Tau) #Iinac,Icout changed def here to instantaneous iloadavg=np.mean(VCout/self.Rload) ploadavg=np.mean(VCout*iloadavg) Ivao=self.get_Ivao(self.get_Vsense(VCout)) #current Ivao vs dt 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 icnt>=0: #self.Clock+=thalf self.addClockThalf() iterdat.append([self.Clock, self.CL0.QCOUT.get_VCout(), self.get_Vss(), IIvao[-1], np.mean(Ivao), self.Vvao, self.Rload, 0.5*ploadavg, np.mean(Pinac), np.mean(Pcout), np.mean(Icout)]) if False: self.gpp["iterLV"].stdin.write("$idat << EOD \n") for i in range(0, len(Tau)): ll=[thalf*Tau[i]/np.pi, 2*Iinac[i], self.get_Vsense(VCout[i]), DVCout[i], 2*Pinac[i], 2*Pcout[i], Duty[i], Vimo[i], Ivao[i], IIvao[i].real] self.gpp["iterLV"].stdin.write(("{} "*len(ll)+" \n").format(*ll)) self.gpp["iterLV"].stdin.write("EOD \n") self.gpp["iterLV"].stdin.write("set multiplot \n") self.gpp["iterLV"].stdin.write("set origin 0,0 \n set size 1,0.5 \n") self.gpp["iterLV"].stdin.write("plot $idat u 1:9 w l lw 2 title 'Ivao' \n") self.gpp["iterLV"].stdin.write("set origin 0,0.5 \n set size 1,0.5 \n") self.gpp["iterLV"].stdin.write("plot $idat u 1:10 w l lw 2 title 'IIvao' \n") self.gpp["iterLV"].stdin.write("unset multiplot \n") #self.gpp["iterLV"].stdin.write("plot $idat u 1:4 w l lw 2 title 'DVCout' \n") #self.gpp["iterLV"].stdin.write("plot $idat u 1:3 w l lw 2 title 'Vsense' \n") #self.gpp["iterLV"].stdin.write("plot $idat u 1:6 w l lw 2, $idat u 1:7 w l lw 2, $idat u 1:8 w l lw 2 \n") gpp=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE) gpp.stdin.write("set term qt size 1900,950 \n") gpp.stdin.write("$idat << EOD \n") for irow in iterdat: gpp.stdin.write(("{} "*len(irow)+" \n").format(*irow)) gpp.stdin.write("EOD \n") gpp.stdin.write("plot $idat u 1:2 w l lw 2 title 'VCout', $idat u 1:8 w l lw 2 title 'ploadavg', $idat u 1:9 w l lw 2 title 'pinacavg', $idat u 1:10 w l lw 2 title 'pcoutavg', $idat u 1:(100*$6) w l lw 2 title '100*Vvao', $idat u 1:(100*$3) w l lw 2 title '100*Vss' \n pause -1 \n") return None ######################################################################################################### # 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 self.Vreg_currsense/self.get_IMO_max() def get_Kdivac(self): return self.RBac/(self.RAac+self.RBac) def get_Vvinac(self, vinac): return vinac*self.get_Kdivac() def get_Vinac_Correction_Factor(self): return 1.0 #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 #ktst=[ [3.0*k/1000, self._get_Kvff(3.0*k/1000)] for k in range(1,1000) ] #for k in ktst: # print(k[0], k[1], k[0]**2/2, (k[0]**2/2)/k[1]) #quit() 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] \ ) #def getLevel(k): # return min(max(1, round( 8-log10(k/2.82)/log10(.856) )), 8) # 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 self.Lx*self.get_Kdivvs()/((0.1e-9)*self.Rsequiv) #appears 100pF cap inside 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): return sqrt(2)*poutLx / ( (self.get_ILxTop(vin, poutLx)**2)*self.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_Rzc_Czc(self): rzc=1/(self.Gmc*abs(self.get_Gpsc(self.Fzc))) czc=1/(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) def addClockThalf(self): self.cntThalf+=1 self.Clock+=1.0/(self.Fac*2.0) 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)) 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_ILx_LoHi_Range(self, nlo, nhi, pout): #NO GOOD FOR DUTY < .5 tpHalf=(1/self.Fpwm)/2 iCpout = -2*pout/self.Vout qCap=self.Cout*self.Vout L0Lst, mLst0, L1Lst, mLst1, InCapLst = [],[],[],[],[] ph1_on, ph1_off, ph1_end, v1, i1, d1, dt1_on, dt1_off, iL1_on, iL1_off = self._get_ILx_Vals(nlo-(1/2), pout) for n in range(nlo, nhi): #get new ph0 ph0_on, ph0_off, ph0_end, v0, i0, d0, dt0_on, dt0_off, iL0_on, iL0_off = self._get_ILx_Vals(n, pout) iL1_at_ph0_on = iL1_on+(iL1_off-iL1_on)*(ph0_on-ph1_on)/(ph1_off-ph1_on) #Iinit + Islope*(t-tinit) iL0_at_ph1_off = iL0_on+(iL0_off-iL0_on)*(ph1_off-ph0_on)/(ph0_off-ph0_on) iL0_at_ph1_end = iL0_on+(iL0_off-iL0_on)*(ph1_end-ph0_on)/(ph0_off-ph0_on) L0Lst.append([ph0_on, iL0_on]) L1Lst.append([ph1_off, iL1_off]) InCapLst.append([ph0_on, iL0_on+iL1_at_ph0_on, iCpout, qCap/self.Cout]); qCap+=(dt1_on-tpHalf)*(iCpout) InCapLst.append([ph1_off, iL1_off+iL0_at_ph1_off, iCpout, qCap/self.Cout]); #InCapLst.append([ph1_off, iL0_at_ph1_off, iCpout+iL1_off, qCap/self.Cout]); qCap+=(dt1_off)*(iCpout+(iL1_off+iL1_on)/2) InCapLst.append([ph1_off, iL1_off+iL0_at_ph1_off, iCpout+iL1_off, qCap/self.Cout]); qCap+=(dt1_off)*(iCpout+(iL1_off+iL1_on)/2) #InCapLst.append([ph1_end, iL0_at_ph1_end, iCpout+iL1_on, qCap/self.Cout]); InCapLst.append([ph1_end, iL1_on+iL0_at_ph1_end, iCpout+iL1_on, qCap/self.Cout]); #add to mlst1 mLst1.append([(ph1_on+ph1_end)/2, v1, i1, d1, dt1_on, dt1_off, iL1_on, iL1_off]) #get new ph1 ph1_on, ph1_off, ph1_end, v1, i1, d1, dt1_on, dt1_off, iL1_on, iL1_off = self._get_ILx_Vals(n+(1/2), pout) iL0_at_ph1_on = iL0_on+(iL0_off-iL0_on)*(ph1_on-ph0_on)/(ph0_off-ph0_on) iL1_at_ph0_off = iL1_on+(iL1_off-iL1_on)*(ph0_off-ph1_on)/(ph1_off-ph1_on) iL1_at_ph0_end = iL1_on+(iL1_off-iL1_on)*(ph0_end-ph1_on)/(ph1_off-ph1_on) L0Lst.append([ph0_off, iL0_off]) L1Lst.append([ph1_on, iL1_on]) InCapLst.append([ph1_on, iL1_on+iL0_at_ph1_on, iCpout, qCap/self.Cout]); qCap+=(dt0_on-tpHalf)*(iCpout) InCapLst.append([ph0_off, iL0_off+iL1_at_ph0_off, iCpout, qCap/self.Cout]); #InCapLst.append([ph0_off, iL1_at_ph0_off, iCpout+iL0_off, qCap/self.Cout]); qCap+=(dt0_off)*(iCpout+(iL0_off+iL0_on)/2) InCapLst.append([ph0_off, iL0_off+iL1_at_ph0_off, iCpout+iL0_off, qCap/self.Cout]); qCap+=(dt0_off)*(iCpout+(iL0_off+iL0_on)/2) #InCapLst.append([ph0_end, iL1_at_ph0_end, iCpout+iL0_on, qCap/self.Cout]); InCapLst.append([ph0_end, iL0_on+iL1_at_ph0_end, iCpout+iL0_on, qCap/self.Cout]); #add to mlst0 mLst0.append([(ph0_on+ph0_end)/2, v0, i0, d0, dt0_on, dt0_off, iL0_on, iL0_off]) result = [ L0Lst, mLst0, L1Lst, mLst1, InCapLst ] #print(result) names = [ "IL0peaks.dat", "IL0vals.dat","IL1peaks.dat", "IL1vals.dat", "IInICapVCap.dat" ] # save them to file here dammit for idx in range(0,len(result)): with open(names[idx],"w") as fout: for ll in result[idx]: for nn in ll: print(nn, " ", end='', file=fout) print("", file=fout) ######################################################################################################### 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 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.LoopDatC=[] self.fcycIvao=int(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.lcDatTau=np.array([]).reshape(5,0) 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.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 = False, 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_IIcaox_Last(self): return self.IIcaox[-1].real def update_Vcaox_with_IIcaox(self): self.Vcaox_prev=self.get_Vcaox() self.Vcaox=max(self.parentCL.Vcaox_bounds[0], min(self.Vcaox + self.get_IIcaox_Last(), self.parentCL.Vcaox_bounds[1])) def get_Vcaox(self): return self.Vcaox def get_PWMx(self): return self.parentCL.Vramp-(self.get_Vcaox() + self.IIcaox.real) #def set_Arr_All(self, ilxarr, flagSavePrev): #def set_Arr_All(self, flagSavePrev): def set_Arr_All(self, nph, tauroot, flagSavePrev): #self.Vcsx=self.get_Vcsx(ilxarr) 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): if self.tStamp != self.parentCL.Clock: print("ERROR: CAx_CSx_Obj arrays not current, call set_Arr_All first. SHOULD NOT HAPPEN DAMMIT") quit() if self.has_Root_PWMx(): return findRoot( lambda t: np.interp(t, self.parentCL.Tau, self.PWMx), tlo, thi ) else: return None def get_FnRoot_PWMx(self, tauroot): return np.interp(tauroot, self.parentCL.Tau, self.PWMx) ######################################################################################################### class ILx_Tau_Obj: def __init__(self, parentcl, qcobj, ilx=0.0): self.parentCL=parentcl 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.parentCL.parent.Lx) self.ILxCurrArr=None def configure(self, nph0, tauOn): #def configure(self, nphx, tauOn): #self.nPh0 = (npwm%self.parentCL.npwmMax)*self.parentCL.ntonPh self.nPh0 = nph0 #self.Clock=self.parentCL.clk self.TauOn=tauOn self.QCobj.configure(self.TauOn) self.ILxCurrArr=None #must regenerate after config self.set_ILxCurrArr(self.parentCL.Tau) 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.parentCL.parent.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.parentCL.parent.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.parentCL.parent.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, flagNegOk=False, 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] ) if flagNegOk: return result else: 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.parentCL.parent.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.parentCL.parent.Lx) 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.parentCL.parent.Lx)*dt + self.Kac*dSin/self.Wac ) def get_DE_dI(self, iI, iF): return 0.5*self.parentCL.parent.Lx*(iF**2 - iI**2) def set_ILx(self, ilx): self.ILx=max(0.0, ilx) def set_ILxBeg_to_ILxEnd(self): self.parentCL.parent.PM0.add_E_T("DE_Vinac_On", self.get_DE_Vinac(0.0, self.TauOn), self.parentCL.tpwm*(self.TauOn)) #self.parentCL.parent.PM0.add_E_T("DE_Vinac_Off", self.get_DE_Vinac(self.TauOn, 1.0, self.parentCL.QCOUT._get_QCout_On(self.TauOn)/self.parentCL.QCOUT.Cout), self.parentCL.tpwm*(1.0-self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_Vinac_Off", self.get_DE_Vinac(self.TauOn, 1.0, self.QCobj.get_VCout_Tau(self.TauOn)), self.parentCL.tpwm*(1.0-self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_Vinac_Tot", self.get_DE_Vinac(0.0, self.TauOn), self.parentCL.tpwm*(self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_Vinac_Tot", self.get_DE_Vinac(self.TauOn, 1.0, self.QCobj.get_VCout_Tau(self.TauOn)), self.parentCL.tpwm*(1.0-self.TauOn)) de_ilx_vcout = self.get_DE_VCout(self.TauOn, 1.0, self.QCobj.get_VCout_Tau(self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", self.get_DE_Vinac(0.0, self.TauOn), self.parentCL.tpwm*(self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", self.get_DE_Vinac(self.TauOn, 1.0, self.QCobj.get_VCout_Tau(self.TauOn)), self.parentCL.tpwm*(1.0-self.TauOn)) self.parentCL.parent.PM0.add_E_T("DE_ILx_Net", de_ilx_vcout, self.parentCL.tpwm*(1.0-self.TauOn)*0) iEnd =self.get_ILx_Tau(1.0) self.set_ILx(iEnd) 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, ilxpk, clk, tauOn): #def configure(self, clk, tauOn): def configure(self, tauOn): #self.ILxPeak=ilxpk #self.Clock=clk self.TauOn=tauOn 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): #must add inductor currents externally, out of phase, etc. return self._get_ICout_On(0*tau) #MESSY 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 = -2*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=(self._get_QCout_On(self.TauOn)**2 + 2.0*de_cout_net*self.Cout) if qnew2>0.0: return np.sqrt(qnew2) else: return 0.0 def get_QCout_Tau(self, tau, flagNegOk=False): #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] ) if flagNegOk: return result else: return np.fmax(result, 0.0) def get_PCout_Lx(self): 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 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 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 set_QCoutBeg_to_QCoutEnd(self): de_ilx_vcout = -2*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.get_QCout_Tau(1.0) ) def get_ICoutCurrArr(self): return self.get_ICout_Tau(self.parentCL.Tau) #----------------------------------------------------------------------------------------------------------------------------------- def stepPWM(self, n0, n, plotdat=False, plotDatTau=False, debug=False, dtau=1.0/(2.0*nptsIcaox), eps=1.0e-6): self.lcDatTau_str= "self.Clock+self.tpwm*self.Tau, self.IL0.ILxCurrArr, self.IL0.get_ILxCout_Arr(), self.CA0.Icaox, self.CA0.IIcaox.real" self.lcDatTau_lblstr=["Tau:ILxCurrArr, Tau:Icaox","Tau:get_ILxCout_Arr(), Tau:IIcaox.real"] 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("lcDatRT") 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, 10*self.CA0.TauRootFin, self.CA0.ErrAvg_PWMx, "#+ \ #"self.StatsCL.nPh0[-1], self.StatsCL.TauOn[-1], self.StatsCL.Vimo[-1], self.StatsCL.ILx_Vinac[-1], self.StatsCL.ILx_VCout[-1], self.StatsCL.VCoutAvg[-1] " 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"] #self.parent.PM0.reset() for ncnt in range(n0, n0+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.Vvao)) #assuming Vvao fluctuations small compared to Vvinac over tpwm if self.CA0.get_Vcaox()>self.Vramp_min: #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 AN INSIDiOUS WAY, not impressed, python FAHHHHKKK self.CA0.solnFound = True self.CA0.set_Arr_All(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.get_TauRoot_PWMx() if self.CA0.has_Root_PeakLimit(): taupklim=self.CA0.get_TauRoot_PeakLimit(0.0, self.CA0.TauRoot) if taupklim is not None and taupklim 0.0): self.CA0.set_Arr_All(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?? 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()) self.CA0.ErrAvg_PWMx=(self.CA0.ErrAvg_PWMx + pwmx_root)/2.0 self.CA0.solnFound = True if not self.CA0.solnFound: if debug or 1: if self.CA0.TauRoot > self.parent.Dmax: print(" OOB: Root exceeds bounds: Dmax limiting, found tau_root=", ncnt, self.CA0.TauRoot, end=" \n") else: print(" No Root found in PWMx: assuming Dmax limited, tau_root=", ncnt, self.CA0.TauRoot, end=" \n") self.CA0.TauRoot=self.parent.Dmax self.CA0.solnFound = True if self.CA0.solnFound: if self.CA0.TauRootFin is None: self.CA0.TauRootFin=self.CA0.TauRoot self.CA0.set_Arr_All(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() self.IL0.set_ILxBeg_to_ILxEnd() #also updates power meter self.QCOUT.set_QCoutBeg_to_QCoutEnd() #also updates power meter self.addClockTPWM() if plotdat>1 or 0: try: self.lcDatRT=list( eval("zip(*["+self.lcDatRT_str+"])") ) self.gppPlotter.gpp_Plot(self.lcDatRT, self.lcDatRT_str, self.lcDatRT_lblstr, self.gppPlotter.gpp["lcDatRT"]) except: pass if plotDatTau: self.lcDatTau[3]=self.lcDatTau[1]+np.interp(self.lcDatTau[0]-self.tpwm/2.0, self.lcDatTau[0], self.lcDatTau[1]) self.lcDatTau[4]=self.lcDatTau[2]+np.interp(self.lcDatTau[0]-self.tpwm/2.0, self.lcDatTau[0], self.lcDatTau[2]) self.gppPlotter.gpp_Plot(self.lcDatTau.transpose(), self.lcDatTau_str, self.lcDatTau_lblstr) if plotdat: self.gppPlotter.gpp_Plot(self.lcDatPks, self.lcDatPks_str, self.lcDatPks_lblstr) 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_ZType2_VAO(self, bk, f): return self.get_Zcomp_type2(bk, f, self.Rzv, self.Czv, self.Cpv) def get_ZType2_VAO_DC(self): #hack that returns approx integration constant for DC component: 1/(cp+cz) return self.get_Zcomp_type2(1, 1/(2*pi*1j), 0, self.Czv, self.Cpv) 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() for k in pfc0.__dict__.keys(): print(k,":",pfc0.__dict__[k]) print("ICout at Pmax:", pfc0.get_ICout(sqrt(2)*pfc0.Vac_rms, 2*pfc0.Ptot), "Vripple:", pfc0.get_Vripple()) print("Vvinac max:", pfc0.get_Vvinac(sqrt(2)*pfc0.Vac_rms), "kvff:", pfc0._get_Kvff(pfc0.get_Vvinac(sqrt(2)*pfc0.Vac_rms))) print("IMO_max:", pfc0.get_IMO_max(), "Vimo_max:", pfc0.Rimo*pfc0.get_IMO_max()) print("Duty at top:", pfc0.get_Duty(sqrt(2)*pfc0.Vac_rms), " at peak voltage:", sqrt(2)*pfc0.Vac_rms) print("ILxTop:", pfc0.get_ILxTop(pfc0.Vac_rms,pfc0.Ptot/2), "dILxTop:", pfc0.get_dILxTop(pfc0.Vac_rms), "ILxPeak:", pfc0.get_ILxPeak()) print("dILxTop/ILxTop:", pfc0.get_dILxTop(pfc0.Vac_rms)/pfc0.get_ILxTop(pfc0.Vac_rms,pfc0.Ptot/2)) print("Rsequiv:", pfc0.Rsequiv, " --> need Vreg_currpeak >= ", pfc0.Vreg_currsense*pfc0.get_ILxPeak()/pfc0.get_ILxTop(pfc0.Vac_rms, pfc0.Ptot/2)) print("LxMin_dILx:", pfc0.get_LxMin_dILx(pfc0.Vac_rms, pfc0.Ptot/2, pfc0.Fpwm), "for fracILx:", pfc0.fracILx) print("LxMin_CCM at Pmax/4:", pfc0.get_LxMin_CCM(pfc0.Vac_rms, pfc0.Ptot/(2*4), pfc0.Fpwm)) print("LxMin_CCM at Pmax/3:", pfc0.get_LxMin_CCM(pfc0.Vac_rms, pfc0.Ptot/(2*3), pfc0.Fpwm)) print("PoutMin_CCM:", pfc0.get_PoutMin_CCM(pfc0.Vac_rms, pfc0.Fpwm), "with Lx:", pfc0.Lx) #print(" phase disable-->1/2, skip half cycle-->1/2; must have PoutMin << PoutMax/2") print("_get_ILx_Vals at n=833:", pfc0._get_ILx_Vals(833, pfc0.Ptot/2.0)) pfc0.Rload=225.0 #pfc0.QCout=340*pfc0.Cout pfc0.CL0.QCOUT.QCout=340*pfc0.CL0.QCOUT.Cout pfc0.Vvao=4.8/3 pfc0.CL0.CA0.Vcaox=4.6/3 pfc0.CL0.IL0.ILx=0.0 iIn, iFin = 0, round(1667) if False or 0: pfc0.CL0.stepPWM(iIn, iFin-iIn, plotdat=1, plotDatTau=1, debug=False) quit() if False or 1: pfc0.QCss=1.8*pfc0.Css pfc0.iter_Loop_Voltage(60, useCurrentLoop=True) quit() pfc0.get_ILx_LoHi_Range(830, 833, 800) gpp=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE) gpp.stdin.write("call 'gpILxPeaks.gp' 'IInICapVCap.dat' 'IL0peaks.dat' 'IL1peaks.dat' 'IL0vals.dat' 'IL1vals.dat' \n") gpp.stdin.flush() ftmp=[p/100.0 for p in range(0, 300, 1)]; flo=1 gptmp=[ pfc0.get_Z_lst(flo*pow(10,f), pfc0.get_GTOTv(flo*pow(10,f))) for f in ftmp] with open("gvtotBode.dat","w") as fout: for t in gptmp: print(t[0], t[2], t[3], file=fout) gpp=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE) gpp.stdin.write("call 'gpBode.gp' 'gvtotBode.dat' \n") gpp.stdin.flush() ftmp=[p/100.0 for p in range(301, 600, 1)] gptmp=[ pfc0.get_Z_lst(pow(10,f), pfc0.get_GTOTc(pow(10,f))) for f in ftmp] with open("gctotBode.dat","w") as fout: for t in gptmp: print(t[0], t[2], t[3], file=fout) gpp=subprocess.Popen(['gnuplot', '-p'], shell=False, universal_newlines=True, stdin=subprocess.PIPE) gpp.stdin.write("call 'gpBode.gp' 'gctotBode.dat' \n") gpp.stdin.flush() ################################## ################################## ################################## 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 #-----------------------------------------------------------------------------------------------------------------------------------------------