連續函數離散化-以SOGI為例

0. 引言

0.1 本文內容

基於SOGI函數，將s域傳遞函數轉換為離散的z域函數，並以m語言形式進行實現，在simulink中封裝為m-function並進行驗證

0.2 學到什么

離散化方法
函數程序實現方法

1. SOGI簡介

以TI官方文檔中單相鎖相環中SOGI應用為例

框圖如下所示

正弦信號經過SOGI可得到同相信號及正交信號

2. 傳遞函數

同相傳遞函數

\[H_{d}(s)=\frac{v^{\prime}}{v}(s)=\frac{k \omega_{n} s}{s^{2}+k \omega_{n} s+\omega_{n}^{2}} \tag{1} \]

正交信號傳遞函數為

\[H_{q}(s)=\frac{q v^{\prime}}{v}(s)=\frac{k \omega_{n} ^2}{s^{2}+k \omega_{n} s+\omega_{n}^{2}} \tag{2} \]

3. 離散化

采用雙線性變換將s域函數離散至Z域

3.1 手動離散

雙線性變換公式為

\[s=\frac{2}{T_{s}}\frac{z-1}{z+1} \tag{3} \]

將式3代入式1得到

\[H_{d}(z)=\frac{k \omega_{n} \frac{2}{T_{s}}\frac{z-1}{z+1}}{(\frac{2}{T_{s}}\frac{z-1}{z+1})^{2}+k \omega_{n} (\frac{2}{T_{s}}\frac{z-1}{z+1})+\omega_{n}^{2}} \tag{4} \]

這里使用以下兩個替換

\[x=2k\omega_{n}T_{s} \tag{5} \]

\[y=( \omega_{n} T_{s})^2 \tag{6} \]

得到

\[H_{d}(z)=\frac{\left(\frac{x}{x+y+4}\right)+\left(\frac{-x}{x+y+4}\right) z^{-2}}{1-\left(\frac{2(4-y)}{x+y+4}\right) z^{-1}-\left(\frac{x-y-4}{x+y+4}\right) z^{-2}}=\frac{b_{0}+b_{2} z^{-2}}{1-a_{1} z^{-1}-a_{2} z^{-2}} \tag{7} \]

同理得到正交函數的離散形式

\[H_{q}(z)=\frac{\left(\frac{k \cdot y}{x+y+4}\right)+2\left(\frac{k \cdot y}{x+y+4}\right) z^{-1}+\left(\frac{k \cdot y}{x+y+4}\right) z^{-2}}{1-\left(\frac{2(4-y)}{x+y+4}\right) z^{-1}-\left(\frac{x-y-4}{x+y+4}\right) z^{-2}}=\frac{q b_{0}+q b_{1} z^{-1}+q b_{2} z^{-2}}{1-a_{1} z^{-1}-a_{2} z^{-2}} \tag{8} \]

3.2 基於MATLAB的離散方法

看完上面的離散過程，很明顯，太麻煩，有沒有簡單點的方法呢？哎，還真有，MATLAB只需要一條命令就能搞定
MATLAB中c2d命令可通過多種離散方法將連續函數離散化，這里為保持一致，同樣以雙線性變換（tustin）為例進行介紹
(了解更多c2d命令，請點擊了解詳情)
具體用法如下

sysd = c2d(sys,Ts,'method')

其中，sys與sysd分別為離散前后函數，Ts為采樣周期，method為離散化方式，這里就是tustin

直接給出離散過程的MATLAB代碼

%%定義s為傳遞函數
s = tf('s');    

%%定義各參數
k = 0.5;
Wn = 100*pi;    %%50Hz
Ts = 1e-4;          %%10kHz

%%寫出傳遞函數
Hd_s = k*Wn*s/(s^2+k*Wn*s+Wn^2);
Hq_s = k*Wn^2/(s^2+k*Wn*s+Wn^2);

Hd_z = c2d(Hd_s,Ts,'tustin')
Hq_z = c2d(Hq_s,Ts,'tustin')

運行結果為

Hd_z =
 
  0.007791 z^2 - 0.007791
  -----------------------
  z^2 - 1.983 z + 0.9844
 
Sample time: 0.0001 seconds
Discrete-time transfer function.


Hq_z =
 
  0.0001224 z^2 + 0.0002448 z + 0.0001224
  ---------------------------------------
          z^2 - 1.983 z + 0.9844
 
Sample time: 0.0001 seconds
Discrete-time transfer function.

3.3 對比

上面已經給出了采用MATLAB進行離散的結果，采用同樣的參數，這里基於式5-8，給出傳統計算方式的結果

Parameter	value	Parameter	value
b0	0.0078	qb0	0.00012238
b1	0	qb1	0.00024476
b2	-0.0078	qb2	0.00012238
a1	1.9834	a2	-0.9844

可能會看到，這里系數正負號與MATLAB計算出結果有所不同，這里實際結果沒錯哈，認為錯了的自己好好檢查！

4.SOGI的程序實現

既然已經得到離散的SOGI函數，如何將其寫成程序呢，這里以MATLAB語言為例，C語言同理

4.1 離散序列的獲得

根據式7和8，我們知道

\[\frac{U_{o}(z)}{U_{i}(z)}=\frac{b_{0}+b_{2} z^{-2}}{1-a_{1} z^{-1}-a_{2} z^{-2}} \tag{9} \]

\[\frac{U_{qo}(z)}{U_{i}(z)}==\frac{q b_{0}+q b_{1} z^{-1}+q b_{2} z^{-2}}{1-a_{1} z^{-1}-a_{2} z^{-2}} \tag{10} \]

容易寫成序列方程

\[\ U_{o} (k)-a_{1}U_{o} (k-1)-a_{2}U_{o} (k-2)=b_{0}U_{i}(k)+b_{2}U_{i}(k-2) \tag{11} \]

\[\ U_{qo} (k)-a_{1}U_{qo} (k-1)-a_{2}U_{qo} (k-2)=qb_{0}U_{i}(k)+qb_{1}U_{i}(k-1)+qb_{2}U_{i}(k-2) \tag{12} \]

4.2 封裝一個m-function

根據上面的式子我們很容易可以寫出相應的程序，但為了在simulink中驗證程序的正確性，我們在這里把SOGI封裝為一個m-function塊以便使用
不了解Matlab的function塊功能的自行百度

很容易知道，對於一個完整的SOGI函數，有一個輸入端，兩個輸出端。函數中各參數均設定為外部給定
下面直接給出相應程序

%%
%%函數聲明
function [uo,quo]   = Orthogonal_Generator(ui,Ts,w,k)

%%
%%定義各中間變量
persistent x;
persistent y;
persistent temp;
persistent b0;
persistent b2;
persistent a1;
persistent a2;
persistent qb0;
persistent qb1;
persistent qb2;
persistent u0;          %%代表ui(k)
persistent u1;          %%代表ui(k-1)
persistent u2;          %%代表ui(k-2)
persistent osg_u0;  %%代表uo(k)
persistent osg_u1;  %%代表uo(k-1)
persistent osg_u2;  %%代表uo(k-2)
persistent osg_qu0; %%代表uqo(k)
persistent osg_qu1; %%代表uqo(k-1)
persistent osg_qu2; %%代表uqo(k-2)

%%
%%初始化各中間變量
if isempty(x)     x= 0; 
end
if isempty(y)     y= 0; 
end
if isempty(temp)  temp= 0; 
end
if isempty(b0)    b0= 0; 
end
if isempty(b2)    b2= 0; 
end
if isempty(a1)    a1= 0; 
end
if isempty(a2)    a2= 0; 
end
if isempty(qb0)   qb0= 0; 
end
if isempty(qb1)   qb1= 0; 
end
if isempty(qb2)   qb2= 0; 
end
if isempty(u0)    u0= 0; 
end
if isempty(u1)    u1= 0; 
end
if isempty(u2)    u2= 0; 
end
if isempty(osg_u0)   osg_u0= 0; 
end
if isempty(osg_u1)   osg_u1= 0; 
end
if isempty(osg_u2)   osg_u2= 0; 
end
if isempty(osg_qu0)   osg_qu0= 0; 
end
if isempty(osg_qu1)   osg_qu1= 0; 
end
if isempty(osg_qu2)   osg_qu2= 0; 
end

%%
%%各系數賦值
x      = 2*k*w*Ts;
y      = w*Ts*w*Ts;
temp   = 1/(x+y+4.0);
b0     = x*temp;
b2     = (-1.0)*b0;
a1     = (2.0)*(4.0-y)*temp;
a2     = (x-y-4)*temp;
qb0    = (k*y)*temp;
qb1    = qb0*(2.0);
qb2    = qb0;

%%
%%計算過程，對應式11離散序列
u0            = ui;
osg_u0        = (b0*(u0-u2)) + (a1*osg_u1) + (a2*osg_u2);
osg_u2        = osg_u1;
osg_u1        = osg_u0;
%%對應式12離散序列
osg_qu0       = (qb0*u0) + (qb1*u1) + (qb2*u2) + (a1*osg_qu1) + (a2*osg_qu2);
osg_qu2       = osg_qu1;
osg_qu1       = osg_qu0;
%%更新序列值
u2            = u1;
u1            = u0;
%%輸出
uo            =osg_u0;
quo           =osg_qu0;

程序有了，我們在simulink中的Library中找到MATLAB Function，寫入上面函數即可

為了進行測試，我們給定一個幅值100，頻率50Hz的正弦信號，其余與上文相同，整個測試模型如下圖所示

同時，要想模型按離散進行仿真，還需要進行相應設置如下圖所示，關鍵在於固定步長

至此，程序編寫及模型搭建，環境搭建就已經完成

4.3 測試

這里運行simulink仿真，將輸入信號，輸出同相信號與輸出正交信號進行對比，如下圖所示

很顯然，在經過兩個周期后，同相輸出信號與輸入重疊，正交信號相差為90°，測試結果表明程序及模型的正確性