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

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

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

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

上述方法均無法滿足實際應(yīng)用對開關(guān)角度求解的實時性、準(zhǔn)確性和可靠性的要求。工業(yè)應(yīng)用中常采用查表法[15]或人工神經(jīng)網(wǎng)絡(luò)法[16-17]來獲取開關(guān)角度。查表法將離線求解的開關(guān)角度按照一定的抽樣間隔存入微控制器中供程序?qū)崟r查詢。此方法需要占用大量芯片上的存儲空間，當(dāng)調(diào)制比位于采樣間隔內(nèi)時，需要通過插值近似確定開關(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)確度，但同時也會增加存儲空間和計算時間。

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

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

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

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

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

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

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

(1)

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

(2)

其中

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

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

(3)

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

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

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

(4)

其中

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

其中

-1

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

已知變元個數(shù)為5的初等對稱多項式定義為

(6)

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

(7)

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

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

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

(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單變元線性方程組中A2，B2，A3，B3，A4，B4，A5，B5是僅與調(diào)制比m相關(guān)的表達式。

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

3 SHE方程組在線求解

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

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

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

(9)

將式(9)展開得

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

(10)

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

(11)

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

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

(12)

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

設(shè)f(x)為實系數(shù)多項式，通過如下方式構(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))為多項式fk-1(x)除fk(x)的余式；Sturm序列的前N項fi(x)≠0(i=0，1，2…，N)為f(x)具有N個實根；fN+1(x)=0為Sturm序列的結(jié)束。

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

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

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

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

圖4 牛頓迭代法與弦截法求根對比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é)所提到的方法進行最終開關(guān)角度的處理：αi=π-βi，i=2，4為下降沿；開關(guān)角為αi，i=1，3，5為上升沿。

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

4.1 案 例

以開關(guān)角個數(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個實根，證明SHE方程組在該調(diào)制比下有解；然后使用二分法求出每個區(qū)間，計算結(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) 使用弦截法精確求解方程的根。對這5個區(qū)間使用弦截法迭代，得出5個具體實根為： -0.962 97，-0.684 6，-0.001 8，0.640 4，0.909 0；最后，將求出的5個根三角反變換，得到開關(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)計算，開關(guān)點數(shù)N=5、6時，SHE方程組在調(diào)制比在(-0.8，0.8)范圍內(nèi)有解；N=7、8時，SHE方程組在調(diào)制比在 (-0.79，0.79)內(nèi)有解。開關(guān)點數(shù)正調(diào)制比表示基波初相位為0，負(fù)調(diào)制比表示基波初相位為π。

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

4.3 實時性分析

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

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

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

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

5 實 驗

為了驗證本文提出的算法的有效性，搭建兩電平逆變器實驗平臺如圖8所示。

圖8 兩電平逆變器實驗平臺Fig.8 Experiment platform of two level inverter

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

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

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

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

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

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

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

6 結(jié) 論

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

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

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