王晨旭，何飛，張奇，潘素娜，楊克虎

1. 中國礦業(yè)大學(xué)(北京) 機(jī)電與信息工程學(xué)院，北京 100083；2. 奧爾堡大學(xué)，丹麥 奧爾堡 9220

采用特定諧波消除脈寬調(diào)制技術(shù)(Selected Harmonic Elimination Pulse Width Modulation，SHEPWM)可以精確地消除逆變器輸出電壓中選定的低次諧波。它具有開關(guān)頻率低、輸出電壓質(zhì)量高等優(yōu)點(diǎn)，被廣泛應(yīng)用于各種大功率場合[1-3]。應(yīng)用SHEPWM技術(shù)需要從一組非線性超越方程組(SHE方程組)中實(shí)時(shí)求解一個(gè)周期的開關(guān)角度，以實(shí)現(xiàn)系統(tǒng)的閉環(huán)控制。由于SHE方程組的非線性特性，開關(guān)角度的求解較為困難，目前常用的求解方法有數(shù)值法[4-6]、智能優(yōu)化算法[7-10]、代數(shù)法[11-14]。

數(shù)值法主要包括牛頓法[4]、多項(xiàng)式同倫算法[5]、沃爾什函數(shù)法[6]等，這類算法具有一定的數(shù)學(xué)理論基礎(chǔ)，易于實(shí)現(xiàn)，求解速度快，但對(duì)初值有依賴性，若初值選取不當(dāng)會(huì)導(dǎo)致迭代時(shí)間長甚至發(fā)散。智能優(yōu)化算法從隨機(jī)解出發(fā)，按照遺傳算法[9]、差分進(jìn)化算法[10]等原理進(jìn)行有限次迭代尋找最優(yōu)解，這類算法對(duì)初值要求低、實(shí)現(xiàn)簡單，但其數(shù)學(xué)基礎(chǔ)薄弱，缺乏普遍意義上的理論分析，容易陷入局部最優(yōu)解，無法保證其收斂性。代數(shù)法主要包括吳方法[11]、結(jié)式消元法[12]、Groebner基算法[13]等，這類算法通過數(shù)學(xué)理論對(duì)SHE方程組進(jìn)行化簡、消元，既不需要初值，又可保證求出全局最優(yōu)解。但代數(shù)法的原理較為復(fù)雜，求解時(shí)占用的計(jì)算機(jī)內(nèi)存多，難以在微控制器上實(shí)現(xiàn)。

上述方法均無法滿足實(shí)際應(yīng)用對(duì)開關(guān)角度求解的實(shí)時(shí)性、準(zhǔn)確性和可靠性的要求。工業(yè)應(yīng)用中常采用查表法[15]或人工神經(jīng)網(wǎng)絡(luò)法[16-17]來獲取開關(guān)角度。查表法將離線求解的開關(guān)角度按照一定的抽樣間隔存入微控制器中供程序?qū)崟r(shí)查詢。此方法需要占用大量芯片上的存儲(chǔ)空間，當(dāng)調(diào)制比位于采樣間隔內(nèi)時(shí)，需要通過插值近似確定開關(guān)角度，給控制系統(tǒng)引入一定的穩(wěn)態(tài)誤差。人工神經(jīng)網(wǎng)絡(luò)法使用離線訓(xùn)練的神經(jīng)網(wǎng)絡(luò)模型生成開關(guān)角度，具有較好的響應(yīng)速度。神經(jīng)網(wǎng)絡(luò)畢竟是一種間接生成開關(guān)角度的方法，不能保證結(jié)果的準(zhǔn)確性，雖然增加神經(jīng)網(wǎng)絡(luò)的規(guī)模可以提高準(zhǔn)確度，但同時(shí)也會(huì)增加存儲(chǔ)空間和計(jì)算時(shí)間。

為了克服現(xiàn)存方法的缺點(diǎn)，筆者提出了一種適用于微控制器的SHE方程組實(shí)時(shí)求解方法，該方法分為離線轉(zhuǎn)換和在線求解兩個(gè)階段。離線轉(zhuǎn)換階段使用對(duì)稱多項(xiàng)式和Groebner基理論對(duì)SHE方程組進(jìn)行轉(zhuǎn)換化簡；在線求解階段使用Sturm定理和數(shù)值法求解具體實(shí)根。此方法結(jié)合了代數(shù)法和數(shù)值法的優(yōu)點(diǎn)，實(shí)現(xiàn)精確、高效的求解，具有以下優(yōu)點(diǎn)：無須給定初值；占用計(jì)算資源少，可在微控制器中實(shí)現(xiàn)；求解效率高，滿足控制系統(tǒng)實(shí)時(shí)性的要求；計(jì)算結(jié)果精確度高?；诖藢?shí)時(shí)求解算法，本文搭建了兩電平逆變器實(shí)驗(yàn)平臺(tái)，驗(yàn)證了此算法的有效性。

1 SHE數(shù)學(xué)模型

圖1為兩電平逆變器拓?fù)浣Y(jié)構(gòu)。該拓?fù)溆扇齻€(gè)橋臂組成，分別對(duì)應(yīng)輸出三相交流電壓；每個(gè)橋臂上下2個(gè)開關(guān)管交替導(dǎo)通，以實(shí)現(xiàn)交流電壓輸出；Udc為直流電源，R為負(fù)載電阻。

圖1 兩電平逆變器拓?fù)浣Y(jié)構(gòu)Fig.1 Topology of two-level inverters

兩電平逆變器的輸出相電壓波形如圖2所示。為了簡化模型，通常假定PWM波形關(guān)于1/4周期偶對(duì)稱、半周期奇對(duì)稱。

圖2 兩電平逆變器SHEPWM輸出電壓波形Fig.2 Output voltage waveform of two-level inverters

對(duì)圖2中的波形進(jìn)行傅里葉級(jí)數(shù)展開，可得各次諧波幅值與開關(guān)角度的關(guān)系：

(1)

SHEPWM技術(shù)的基本思想是通過控制開關(guān)角度αi，輸出期望的基波幅值，并使得某些特定次諧波的幅值等于零。以開關(guān)角N=5為例，在控制基波幅值的同時(shí)消除第3、5、7、9次諧波，根據(jù)式(1)得SHE方程組：

(2)

其中

式中，m為調(diào)制比。

調(diào)制比與基波幅值U的關(guān)系定義為

(3)

2 SHE方程組的離線轉(zhuǎn)換

由式(2)可知，SHE方程組為非線性超越方程組，直接求解較為困難，利用對(duì)稱多項(xiàng)式和Groebner基理論，可以將SHE方程組的求解等價(jià)轉(zhuǎn)化為一個(gè)線性方程組和一個(gè)一元高階多項(xiàng)式的求解[19-20]。

首先，對(duì)式(2)中的負(fù)余弦項(xiàng)應(yīng)用變量代換αi=π-βi，使余弦項(xiàng)的符號(hào)統(tǒng)一為正號(hào)：

(4)

其中

采用三角函數(shù)倍角公式將式(4)展開，并作變量代換cosαi=xi，i=1，3，5，cosβi=xi，i=2，4，則SHE方程組被轉(zhuǎn)化為一個(gè)對(duì)稱的多項(xiàng)式系統(tǒng)：

其中

-1

式(5)為一個(gè)9次的多變?cè)匠探M，若使用代數(shù)法直接求解，隨著開關(guān)角個(gè)數(shù)的增加，計(jì)算量和內(nèi)存消耗會(huì)急劇增加。為了降低計(jì)算量，可使用初等對(duì)稱多項(xiàng)式理論將式(5)化簡降階，詳細(xì)的數(shù)學(xué)轉(zhuǎn)換原理和推導(dǎo)過程可參考文獻(xiàn)[13]，這里僅作簡要描述。

已知變?cè)獋€(gè)數(shù)為5的初等對(duì)稱多項(xiàng)式定義為

(6)

使用初等對(duì)稱多項(xiàng)式將式(5)降階，此過程可通過調(diào)用Mathematica軟件中的Symmetric Reduction命令完成。由于式(7)中的表達(dá)式較長，故省略了中間部分：

(7)

f5(x)和f5(e)的階次比較見表1。相比式(5)，式(7)中變量的階次得到了大幅降低，因此計(jì)算量也大大降低。

表1 f5(x)和f5(e)的階次比較

進(jìn)一步計(jì)算式(7)中多項(xiàng)式組的Groebner基[20]，可以得到等價(jià)方程組：

(8)

其中

A2=4m12-48m11+54m10+1 290m9-4 185m8-

9 180m7+47 250m6+14 175m5-198 450m4+28 350m3+340 200m2-42 525m-170 100

B2=36(m10-10m9+255m7-525m6-1575m5+

4 725m4+1 575m3-9 450m2+4 725)

A3=m13-13m12+12m11+462m10-1 515m9-

4 995m8+25 785m7+9 450m6-155 925m5+

99 225m4+283 500m3-255 150m2-127 575m+

127 575

B3=72(m10-10m9+255m7-525m6-1 575m5+

4 725m4+1 575m3-9 450m2+4 725)

A4=m14-14m13+7m12+672m11-2 226m10-

9 870m9+51 975m8+32 130m7-423 360m6+

198 450m5+1 289 925m4-793 800m3-

1 786 050m2+595 350m+893 025

B4=1 008(m10-10m9+255m7-525m6-

1 575m5+4 725m4+1 575m3-9 450m2+

4 725)

A5=m15-15m14+945m12-3 150m11-18 270m10+

97 650m9+99 225m8-1 091 475m7+

396 900m6+5 060 475m5-4 465 125m4-

8 930 250m3+8 930 250m2+4 465 125m-

4 465 125

B5=30 240(m10-10m9+255m7-525m6-

1 575m5+4 725m4+1 575m3-9 450m2+4 725)

可見，e1，e2，…，e5單變?cè)€性方程組中A2，B2，A3，B3，A4，B4，A5，B5是僅與調(diào)制比m相關(guān)的表達(dá)式。

需要指出，上述轉(zhuǎn)化過程的運(yùn)算量較大，需要在計(jì)算機(jī)中借助Maple或者M(jìn)athematica軟件實(shí)現(xiàn)。離線轉(zhuǎn)換結(jié)束后將方程組(8)存入微控制器中，以實(shí)現(xiàn)方程組的在線求解。

3 SHE方程組在線求解

在線求解階段，首先需要根據(jù)給定的調(diào)制比m求解方程組(8)得到e，然后根據(jù)e求解方程組(6)得到x，進(jìn)而得到開關(guān)角度α。其關(guān)鍵和難點(diǎn)在于如何求解方程組(6)。文獻(xiàn)[21]通過結(jié)式消元法將式(6)三角化后再迭代求解，具有轉(zhuǎn)換過程復(fù)雜、中間變量多等缺點(diǎn)。本文利用一元高階多項(xiàng)式根與系數(shù)的關(guān)系，將式(6)的求解轉(zhuǎn)化為一元高階方程的求解，然后利用Sturm定理[22]將一元高階方程的實(shí)根隔離在不同區(qū)間，最后使用弦截法求出精確實(shí)根。

假設(shè)多項(xiàng)式f(x)的實(shí)根是x1，x2，…，xn，則f(x)的形式為

f(x)=(x-x1)(x-x2)…(x-xn)

(9)

將式(9)展開得

f(x)=xn+e1xn-1+…+en-1x+en

(10)

對(duì)比式(9)和式(10)，多項(xiàng)式的系數(shù)與其實(shí)根具有如下關(guān)系：

(11)

通過對(duì)比式(11)和式(6)發(fā)現(xiàn)，在已求得e的前提下，通過構(gòu)造多項(xiàng)式方程有

f(x)=x5-e1x4+e2x3-e3x2+e4x-e5=0

(12)

由式(2)求出x，由多項(xiàng)式方程組(6)的求解轉(zhuǎn)化為式(12)一元高階方程的求解。當(dāng)一元高階方程的階次小于等于4時(shí)，可直接使用求根公式求解；當(dāng)階次大于4時(shí)，可先找到若干互不重疊的區(qū)間，其中每一個(gè)區(qū)間包含一個(gè)實(shí)根，然后對(duì)每一個(gè)區(qū)間使用數(shù)值法求出精確實(shí)根。實(shí)根的隔離通常利用Sturm定理結(jié)合二分法完成，下面簡要介紹Sturm定理。

設(shè)f(x)為實(shí)系數(shù)多項(xiàng)式，通過如下方式構(gòu)造其Sturm序列：

f0(x)=f(x)≠0f1(x)=f′(x)≠0f2(x)=-rem(f0(x)，f1(x))≠0 ?fk+1(x)=-rem(fk-1(x)，fk(x))≠0 ?fN(x)=-rem(fN-2(x)，fN-1(x))≠0fN+1(x)=-rem(fN-1(x)，fN(x))=0

(13)

式中，rem(fk-1(x)，fk(x))為多項(xiàng)式fk-1(x)除fk(x)的余式；Sturm序列的前N項(xiàng)fi(x)≠0(i=0，1，2…，N)為f(x)具有N個(gè)實(shí)根；fN+1(x)=0為Sturm序列的結(jié)束。

Sturm定理：設(shè)(a，b)為任一區(qū)間，若f(a)f(b)≠0，則Sturm序列f0(a)，…，fN(a)的變號(hào)次數(shù)V(a)與Sturm序列f0(b)，f1(b)，…，fN(b)的變號(hào)次數(shù)V(b)的差V(a)-V(b)恰好是f(x)在區(qū)間(a，b)內(nèi)實(shí)根的個(gè)數(shù)。

應(yīng)用Sturm定理，可以判斷多項(xiàng)式在任意區(qū)間內(nèi)的實(shí)根個(gè)數(shù)。由于式(12)的實(shí)根限定在區(qū)間(-1，1)，將此區(qū)間不斷二分，并對(duì)每一個(gè)區(qū)間使用Sturm定理判斷其所包含的實(shí)根個(gè)數(shù)，直到找到N個(gè)互不重疊的區(qū)間且每個(gè)區(qū)間都包含一個(gè)實(shí)根。算法流程如圖3所示。

圖3 Sturm定理應(yīng)用流程Fig.3 The flow chart of using Sturm theorem

根據(jù)Sturm定理得到多項(xiàng)式的所有實(shí)根區(qū)間后，可以用弦截法、牛頓法等[23]數(shù)值法求解方程的實(shí)根。由于每個(gè)實(shí)根的區(qū)間已經(jīng)確定，需要保證求解時(shí)的迭代過程始終位于區(qū)間內(nèi)。本文選用弦截法，而不是牛頓法。因?yàn)榕ｎD法在迭代過程中可能超出區(qū)間范圍，造成求解失敗。如圖4所示，x2是方程在區(qū)間(a1，a2)內(nèi)的一個(gè)實(shí)根，若選x0為初始值，牛頓迭代法將會(huì)錯(cuò)誤地找到區(qū)間外的x3為根，造成求解失敗?？紤]上述情況，選擇弦截法求解，如求解區(qū)間(an-1，an)內(nèi)的實(shí)根時(shí)，弦截法以區(qū)間端點(diǎn)為弦，保證了所有的迭代過程均在區(qū)間內(nèi)完成，非常適合當(dāng)前的求解要求。

圖4 牛頓迭代法與弦截法求根對(duì)比Fig.4 Comparison between the Newton-iteration method and the chord method

通過上述方法求解得到xi后，由三角函數(shù)反變換αi=arccosxi，i=1，3，5，βi=arccosxi，i=2，4得到開關(guān)角度，然后按照第2節(jié)所提到的方法進(jìn)行最終開關(guān)角度的處理：αi=π-βi，i=2，4為下降沿；開關(guān)角為αi，i=1，3，5為上升沿。

4 計(jì)算案例與結(jié)果分析

4.1 案 例

以開關(guān)角個(gè)數(shù)N=5、調(diào)制比m=0.8為例，詳細(xì)說明上述求解方法。

(1)將m=0.8代入式(8)，可得：

(14)

(2)根據(jù)式(12)構(gòu)造一元高次方程：

f(x)=x5+0.1x4-1.311 36x3-0.064 802 6x2+ 0.383 773x+0.000 725 084

(15)

(3)構(gòu)造式(15)的Sturm序列：

(16)

(4)利用Sturm序列劃分解的區(qū)間。把x=-1代入Sturm序列可得V(-1)=5：

(17)

同理，把x=1代入Sturm序列可得V(1)=0：

(18)

因此，方程在(-1，1)區(qū)間內(nèi)含有5個(gè)實(shí)根，證明SHE方程組在該調(diào)制比下有解；然后使用二分法求出每個(gè)區(qū)間，計(jì)算結(jié)果如下：

(-1.000 0，-0.750 0)，(-0.750 0，-0.312 5)，(-0.312 5，0.343 8)，(0.343 8，0.671 9)，(0.671 9，1.000 0)

(5) 使用弦截法精確求解方程的根。對(duì)這5個(gè)區(qū)間使用弦截法迭代，得出5個(gè)具體實(shí)根為： -0.962 97，-0.684 6，-0.001 8，0.640 4，0.909 0；最后，將求出的5個(gè)根三角反變換，得到開關(guān)角度為：15.639°↑，24.626°↓，46.790°↑，50.171°↓，89.892°↑。

4.2 全調(diào)制比內(nèi)的求解結(jié)果

按照文中所述的過程繪制全調(diào)制比內(nèi)的根軌跡圖，如圖5至圖8所示。經(jīng)計(jì)算，開關(guān)點(diǎn)數(shù)N=5、6時(shí)，SHE方程組在調(diào)制比在(-0.8，0.8)范圍內(nèi)有解；N=7、8時(shí)，SHE方程組在調(diào)制比在 (-0.79，0.79)內(nèi)有解。開關(guān)點(diǎn)數(shù)正調(diào)制比表示基波初相位為0，負(fù)調(diào)制比表示基波初相位為π。

圖5 不同開關(guān)點(diǎn)數(shù)時(shí)SHE方程組的解軌跡Fig.5 Roots locus map of SHE equations for different number of switching angles

4.3 實(shí)時(shí)性分析

圖6為在STM32F407控制器上求解SHE方程組所需的時(shí)間曲線。當(dāng)N=5時(shí)，在全調(diào)制比內(nèi)所需的平均時(shí)間為2.25 ms；N=6、7時(shí)，所需的平均時(shí)間分別為3.5 ms和5.25 ms；在N=8時(shí)，所需平均時(shí)間約7 ms，基本滿足控制器實(shí)時(shí)性需求。

圖6 N=5、6、7、8時(shí)該算法的執(zhí)行時(shí)間曲線Fig.6 Executing time curve of proposed method for N=5、6、7、8

圖7為本文提出的實(shí)時(shí)求解算法與牛頓迭代法在STM32F407上的計(jì)算時(shí)間對(duì)比。傳統(tǒng)的數(shù)值算法直接求解SHE方程組，而實(shí)時(shí)求解算法是將SHE方程組轉(zhuǎn)化為一元高階方程求解。由圖7可知，牛頓迭代法所需的時(shí)間約40 ms，實(shí)時(shí)求解法所需時(shí)間約2 ms，求解速度提高了20倍。

圖7 實(shí)時(shí)求解法與牛頓迭代法的計(jì)算時(shí)間對(duì)比Fig.7 Computation time of real-time method and newton-iteration method

5 實(shí) 驗(yàn)

為了驗(yàn)證本文提出的算法的有效性，搭建兩電平逆變器實(shí)驗(yàn)平臺(tái)如圖8所示。

圖8 兩電平逆變器實(shí)驗(yàn)平臺(tái)Fig.8 Experiment platform of two level inverter

控制芯片型號(hào)為ARM Cortex-m 3、處理器STM32F407，開關(guān)器件選用IPM模塊STGIPS30C60，直流側(cè)電壓為30 V，輸出的基波頻率為50 Hz。

實(shí)驗(yàn)選取3組不同開關(guān)角個(gè)數(shù)與調(diào)制比進(jìn)行驗(yàn)證，求解出的開關(guān)角度見表2。

表2 三組開關(guān)角

圖9至圖11為三組開關(guān)角的輸相電壓波形和傅里葉分析。由圖可知，在N=7時(shí)，期望消去的第3、5、7、9、11、13次諧波均已消除；在N=8時(shí)，期望消去的第3、5、7、9、11、13、15次諧波均已消除。說明求解的開關(guān)角度是精確有效的。

圖9 N=7，m=-0.78時(shí)的相電壓 輸出波形和傅里葉分析Fig.9 Output phase voltage waveform and FFT analysis for N=7，m=-0.78

圖10 N=8，m=0.5時(shí)的相電壓輸出波形和傅里葉分析Fig.10 Output phase voltage waveform and FFT analysis for N=8，m=0.5

圖11 N=8，m=-0.74時(shí)的相電壓 輸出波形和傅里葉分析Fig.11 Output phase voltage waveform and FFT analysis for N=8，m=-0.74

6 結(jié) 論

(1) 本文提出了一種適用于微控制器的開關(guān)角實(shí)時(shí)求解方法，在STM32F407處理器中求解5個(gè)開關(guān)角所需時(shí)間為2.25 ms，求解8個(gè)開關(guān)角所需時(shí)間為7 ms，與牛頓迭代法相比，速度提升了20倍。

(2) 該方法很好地解決了傳統(tǒng)數(shù)值算法需要給定初值、智能優(yōu)化算法和代數(shù)法計(jì)算量大等缺陷，具有數(shù)學(xué)理論基礎(chǔ)完備、求解精確度高等優(yōu)點(diǎn)。

(3) 實(shí)驗(yàn)結(jié)果顯示，利用該方法所求出的開關(guān)角能夠精確地控制基波幅值，同時(shí)可以消除指定的各次諧波，驗(yàn)證了算法的正確性與有效性。