霍龍

（沈陽工程學(xué)院 遼寧 沈陽 100136）

二階電路的暫態(tài)分析中，當(dāng)電路處于欠阻尼狀態(tài)時其響應(yīng)為一個衰減振蕩波。通過數(shù)字示波器或數(shù)字仿真可以獲得振蕩波的一組離散化數(shù)據(jù)，它是電路變量的精確時域解。當(dāng)振蕩波是以一組離散的數(shù)據(jù)點給出時，如何準(zhǔn)確地從數(shù)據(jù)中提取振蕩頻率具有重要意義。在電力系統(tǒng)故障分析中，故障電流流波形就是一個衰減振蕩波，在頻域上表現(xiàn)為主頻的諧波形式，利用這些頻率可以進(jìn)行故障測距[1-2]，可以分析暫態(tài)過程對繼電保護及自動裝置工作性能的影響[3]。

以往振蕩頻率測量多采用測量相鄰波峰幅值的方法 由于取用的波形特征數(shù)據(jù)偏少，讀數(shù)誤差較大，計算精度低。并且波形上的其他數(shù)據(jù)點均未取用，造成數(shù)據(jù)的浪費。

為了充分利用衰減振蕩曲線的全部離散化數(shù)據(jù)（數(shù)據(jù)量由采樣間隔決定），文中采用非線性最小二乘優(yōu)化法，通過曲線擬合求出振蕩頻率。首先得出曲線的含有特征參數(shù)的數(shù)學(xué)模型，然后構(gòu)造二乘殘差函數(shù)，設(shè)置特征參數(shù)的初始值，最后通過MATLAB編程用最優(yōu)化算法搜索函數(shù)的極小值。搜索成功后得到擬合曲線的特征參數(shù)，從特征參數(shù)中獲得振蕩頻率的精確值。該算法執(zhí)行速度快，程序通用性強，能夠準(zhǔn)確提取振蕩頻率。

1 原理與算法

二階電路的暫態(tài)過程在欠阻尼時的響應(yīng)為衰減振蕩曲線，其解析式如下：

其中 δ為衰減系數(shù)、ω 為振蕩角頻率、A=f（0）為與 f（t）同量綱的幅值、β為初相角。這4個參數(shù)作為特征參數(shù)唯一地決定了響應(yīng)曲線的形狀（響應(yīng)的表達(dá)式）。當(dāng)響應(yīng)曲線是以離散的數(shù)據(jù)點（tk，fk）給出時，文中擬通過非線性最小二乘優(yōu)化方法確定這4個參數(shù)，進(jìn)而從參數(shù)中得到振蕩頻率

設(shè)擬合曲線模型為：

其中 ai（i=1，2，3，4）為待求參數(shù)。 首先由 n 對離散數(shù)據(jù)點（tk，fk）（k=1，2，…，n）構(gòu)造二乘殘差目標(biāo)函數(shù)：

最小二乘曲線擬合的目標(biāo)就是確定參數(shù)ai，使得目標(biāo)函數(shù)Twoexps的值為最小[4-5]。

在非線性最小二乘估計算法中，有效而完善的算法當(dāng)推Marquadst算法[6]。它是專門為參數(shù)估計的最小二乘準(zhǔn)則而設(shè)計的。Marquadst算法認(rèn)為：1）關(guān)于被估參數(shù)的最小二乘函數(shù)，在接近最小值點時，呈拋物線形態(tài)，且可用Taylor級數(shù)的前三項足夠好地近似，宜取較小步長；2）在遠(yuǎn)離最小值點時，變化陡峭，適于最速下降法處理，步長宜取大。

設(shè)被估參數(shù)向量為a，a的初始近似值為a0，a與初值a0之差為Δa，則a按下式更新：

其中

H為Hessian矩陣，元素的近似值為：

用Δa修正初始值a0得到更新值a。即對于指定精度ε（一般取 0.001），當(dāng) Δa ＞ε 時，用a置換a0重新求解（5）式并代入（4）式，如此反復(fù)更新a值（迭代），直至得到 Δa <ε 時的a值為最終結(jié)果。

算法流程如下：

1）輸入 n 對離散數(shù)據(jù)組（tk，fk），構(gòu)造目標(biāo)函數(shù) Twoexps，設(shè)置精度ε；

2）由數(shù)據(jù)（tk，fk）確定初始值 a0；

3Δ）由式（6）計算 Hessian 矩陣，由（4）式計算 χ2（a）并求得梯度a（χ2）；

4）解矩陣方程式（5）求得修正值Δa；

5）用Δa修正初始值a0得到更新值a；

6）若 Δa ＞ε，用更新值 a 取代 a0重復(fù)步驟 3）、4）、5）直到 Δa <ε，終止更新。

2 初始值的計算

初始值的設(shè)定關(guān)系到最優(yōu)化的收斂速度和結(jié)果的穩(wěn)定性。為避免人工設(shè)置初始值，則要求給出由離散數(shù)據(jù)組（tk，fk）直接確定初始值的方法。 設(shè) a0=[a10，a20，a30，a40]T，本文給出如下計算方法：

參數(shù)a10是A的初始值，可以取yk的最大值，即a10=max（fk），并以此值作為f1，如圖1所示，由此形成模型函數(shù)為衰減振蕩的數(shù)據(jù)組。

圖1 衰減振蕩曲線Fig.1 Damped oscillation characteristic

參數(shù)a30是ω的初始值，設(shè)f（t）的兩個相鄰的過零點分別為 tA和 tB（見圖 1），則：

考慮到離散數(shù)據(jù)點不一定恰好過零，為了求tA和tB的近似值，本文采用直線求交法。以tA為例，取f值異號的相鄰兩個數(shù)據(jù)點 fi-1、 fi（即 fi-1fi＜0），則必有 ti-1＜tA＜ti，連接 fi-1、 fi的直線與時間軸的交點作為tA的近似值。同理，取fj-1、fj的連線與時間軸的交點作為tB的近似值[7]，于是得式（9）：

a40是初相角β的初始值，取決于坐標(biāo)原點的選取。設(shè)f（t）第一零點為 tA，則有 ωtA=π-β，得：

a20是衰減系數(shù)δ的初始值，由于p=δ+jω （也稱固有頻率），有，得：

若離散數(shù)據(jù)（tk，fk）中含有兩個相鄰的極值點（如圖1的fm1和 fm2），則衰減系數(shù) a20也可由式（12）計算[8]：

3 模擬算例和實驗結(jié)果

文中用MATLAB語言編程實現(xiàn)上述算法。下面以算例和實驗說明算法的有效性和準(zhǔn)確性。

例1設(shè)二階電路的暫態(tài)解表達(dá)式為：

為了模擬實測情況，在上述信號中加入1%的白噪聲信號。t在[0，15 ms]中取值，采樣間隔為1 ms，優(yōu)化的終止容差取ε=10-4。

計算結(jié)果如表1所示，角頻率a3的真值為628.32 rad/s，優(yōu)化值為629.039 rad/s。振蕩頻率f的真值為f=100 Hz，優(yōu)化值為：

相對誤差為0.11%，這是從1%的噪聲污染數(shù)據(jù)中優(yōu)化的結(jié)果。若噪聲降為0.5%，則a3的優(yōu)化值為628.678 6，相對誤差為0.057%?？梢姡疚乃惴ㄔ谠紨?shù)據(jù)趨于準(zhǔn)確的前提下，優(yōu)化精度會顯著提高。

表1 例1的計算結(jié)果Tab.1 Computation results of example 1

實驗電路采用Vpp=2 V，f=1 000 Hz的方波作為激勵。為了降低電容中損耗電阻的影響[9]，電容取C=0.01 μF，電感取L=15 mH（損耗電阻rL=17.33 Ω），串聯(lián)可變電阻R可以改變電路的工作狀態(tài)，當(dāng)R=592.53 Ω（含rL）時得到圖2所示波形。

圖2 實驗電壓波形Fig.2 The experimental voltage waveform

數(shù)據(jù)采集時間取 t=（0～110）μs，采樣間隔為 Δt=4 μs。 取ε=10－4時，角頻率a3的優(yōu)化結(jié)果如表2所示，相對誤差為0.000 607 97（約 0.061%）。

表2 例2的計算結(jié)果Tab.2 Computation results of example 2

4 結(jié) 論

針對離散數(shù)據(jù)，本文采用非線性最小二乘優(yōu)化方法進(jìn)行暫態(tài)波形的分析?？梢杂嬎銜簯B(tài)波形中的角頻率、衰減系數(shù)、幅值和初相角，其相對誤差小于10－3。算法易于理解，編程代碼簡明，程序運行可靠（均能進(jìn)入二乘殘差函數(shù)的淺谷中搜索出極小值點）。只要輸入數(shù)據(jù)準(zhǔn)確，采樣區(qū)間?。?.5～2）倍振蕩周期，優(yōu)化結(jié)果可以達(dá)到很高的精度，程序具有很強的實用性和通用性。

[1]王興國，黃少峰.基于復(fù)解析帶通波器的固有頻率自適應(yīng)提取原理和方法[J].電工技術(shù)學(xué)報，2009，24（12）：179－183.WANGXing-guo，HUANGShao-feng.Naturalfrequencyadaptive extracting principle and method based on multiple analysis band-pass filter[J].Transactions of China Electrotechnical Society，2009，24（12）:179－183.

[2]王群，王兆安.時域中非正弦周期電流的分解及其各分量的測量[J].儀器儀表學(xué)報，2000，21（4）：371－875.WANG Qun，WANG Zhao-an.Decomposition in time domain for nonsinusoidal period current and its measurement of each component[J].Chinese Journal of Scientific Instrument，2000，21（4）:371－875.

[3]李庚銀，陳志業(yè)，楊峰.電力系統(tǒng)暫態(tài)波形分析方法[J].中國電機工程學(xué)報，1995，15（3）：204－209.LI Geng-yin，CHEN Zhi-ye，YANG Feng.Analysis method of transient waveform in electric power system[J].Proceedings of the CSEE，1995，15（3）：204－209.

[4]彭景斌，葉進(jìn)寶，王雪嬌.暫態(tài)混沌神經(jīng)網(wǎng)絡(luò)及其在優(yōu)化問題中的應(yīng)用研究[J].現(xiàn)代電子技術(shù)，2009（4）：76－79.PENG Jing-bin,YE Jin-bao,WANG Xue-jiao.Transient chaotic neural network and its optimization of the applied research[J].Modern Electronic Technique,2009（4）：76－79.

[5]陸建明.基于最小二乘虛擬陣元的解模糊方法[J].電子科技，2009（11）：9－11,15.LU Jian-ming.Fuzzy solution based on least squares virtual array algorithm[J].Electronic Science and Technology,2009（11）：9－11,15.

[6]張志涌.精通MATLAB6.5版[M].北京：北京航空航天大學(xué)出版社，2003.

[7]霍龍，俞俊民，于佳.基于曲線求交算法的非線性電阻電路數(shù)值計算[J].哈爾濱理工大學(xué)學(xué)報，2010，15（6）:17－20.HUO Long，YU Jun-min，YU Jia.The Numerical Value Solution of Nonlinear Resistor Circuit Based on Algorithm of Intersection[J].Journal of Harbin University of Science and Technology，2010，15（6）:17－20.

[8]王勤，余定鑫.電路實驗與實踐[M].北京：高等教育出版社，2004.

[9]王璐，楊百瑞.“RLC串聯(lián)電路暫態(tài)過程的研究”實驗中電容系統(tǒng)誤差的測量與修正[J].大學(xué)物理，2008，27（2）:48－49.WANG Lu，YANG Bai-rui.The measurementsofthe systemic uncertainty ofthe capacitorand the related corrections in the “Studies on the transient processes in a RLC series circuit” experiment[J].College Physics，2008，27（2）:48－49.