再來看回授通路:透過採樣電阻來擷取任意兩相電流,根據基爾霍夫電流定律可以算出第三相電流,將三向電流經由Clark變換轉換成 和 ,再透過Park變換轉換成 和 ,作為回授傳入PID控制器構成電流環。
上圖的Park變換和反Park變換需要目前的角度作為輸入;速度PID需要速度作為回饋。所以需要獲得馬達的速度與角度。角度和速度的獲取方法分為有感和無感。有感方式使用霍爾元件(Hall Sensor),安裝在馬達上就可以偵測馬達磁鐵的位置。無感使用觀測器(observer)獲得角度速度訊息,本文將使用擴展卡爾曼濾波觀測器(EKF),輸入為 、 、 、 。使用無感方式不需要霍爾感測器,可以減少連接線的數量,也可以減少成本。
為什麼要使用座標變換?馬達控制大多是在控制速度/轉矩,需要用PID閉環控制正弦交流電壓的振幅和角度,不是很容易實現,所以透過座標變化把正弦交流資訊分解成角度資訊(Q軸控制轉矩)和幅值資訊(D軸控制磁場)單獨控制。
FOC的變換中要滿足等振幅變換,即變換前後振幅不變。
座標變換都分為正向變換和反向變換,正向變換都是對電流進行操作的,反向變換都是對電壓進行操作的。
以下的變換均以聯立和矩陣兩種形式表示,以方便使用。
Clark變換 Clark變換實現了三向座標系 與直角座標系 間的轉換。
7C602099D7A9549EDA9C24D13AB98072
正向Clark變換 Clark反變換則將兩相信號轉換為三相信號。已知三相座標系 ,這三個基底向量不是正交的,所以可以將其正交化為一個直角座標系,命名為 座標系,變換公式為:
表示為矩陣形式:
我們一般在電路中只採集兩相電流,第三相電流可以使用基爾霍夫電流定律來得出( ),故上式也可整理為以下形式:
矩陣形式:
由於變換前後 和 幅值要相同,所以要進行等幅值變換,變換係數為
,即變為
矩陣形式:
(這裡的係數在後文SVPWM里相電壓的振幅與電壓空間向量之間有一個
的係數相抵消)
MATLAB實作為:
image-20220330103933845
反向Clark變換 反Clark變換則將三相信號轉換為兩相信號。根據圖可以寫出:
矩陣形式:
MATLAB實作為:
image-20220330103915885
Park變換 Park變換實現了兩相座標系 與轉子座標系 間的轉換,此變換可以將正弦變數線性化。其中 指向轉子中心, 指向切線方向, 是轉子當前的角度,也就是說 座標系始終跟著轉子同步旋轉。
A20BC4DCAD2986A86187F746C2EDE88A
正向Park變換 則根據上圖可以寫出
很明顯上述變換可以用旋轉矩陣來表示,使用矩陣形式可以很方便地寫出:
(如果 軸為0,則功率全部輸出在 軸上。 )
MATLAB實作為
image-20220330104003922
反Park變換 根據上面的推導可以求反Park變換
同理,使用旋轉矩陣可以求反變換的係數矩陣:
MATLAB實作為
image-20220330103857619
MATLAB仿真 為了更清楚地仿真,這裡不用矩陣形式表示,如需矩陣形式可以看我寫的另一篇文章。
請注意,正向變換都是對電流進行操作的,反向變換都是對電壓進行操作的。但在這節的模擬中,把正變換和反變換連在一起,這樣做沒有實際意義,只是為了驗證變化演算法。
輸入Vd為0,Vq為1,角度為由0到2pi的連續值。
image-20220330103716603
再來看子模組內部,輸入經過兩個逆變換,再經過兩個正變換。
image-20220330103816443
運行查看波形(新版本MATLAB常值輸入為一個圓圈)
image-20220330104412589
SVPWM 根據BLDC的六步換向,可以將一圈分為六個扇區,前文FOC引入章節已經講過,只需要控制每個狀態通電的時間就可以控制轉子到達任意角度。這就是SVPWM。
SVPWM的輸入為 和 ,輸出為三向計數器的比較值。所以應該先判斷用哪兩個相鄰向量,然後計算兩個相鄰向量的作用時長,然後將作用時長轉換成計數器的比較數值送入定時器。以下對這三個步驟進行講解。
扇區判斷 三相電壓可以表示為( 為電壓幅值):
將其轉換為 座標系,可以算出
從這個式子發現,可以從中算出角度資訊從而可以判斷在哪個扇區
由於除法和反三角函數對MCU來說計算量比較大。我們來找一個簡單演算法。
是關於cos的三角函數, 是關於sin的三角函數,可以得到:
扇區 扇區 扇區 扇區 扇區 扇區
可以看到,透過 的正負可以判斷是1/2/3扇區還是4/5/6扇區。
這個條件只將扇區分為兩個部分,我們還需要幾個條件來更細緻地分。將每個磁區的反三角函數範圍計算出來:
扇區
扇區
扇區
扇區
扇區
扇區
觀察這個結論, 和 似乎有關係,回顧反Clark變換, 和 的式子就是這種關係。所以可以把上面的結論往反Clark變換上湊。看一下反Clark變換的圖像,注意我們需要的關係裡的 是乘在 上的,所以我們把 和 反一下,對應的公式為:
產生上述公式的影像,其中黃色的線為 為 ,藍色的線 為
,紅色的線 為
,
image-20220330145027064
可以看到在每個磁區總有一向大於0,兩向小於0,所以 和 的正負可以當作判斷條件之一。我們順便還又一次得到了 這個判斷條件。整理一下上面的式子
扇區
扇區
扇區
扇區
扇區
扇區
所以透過計算
和
的正負可以判斷出是1/6扇區,3/4扇區,2扇區,5扇區的哪一組。
綜上,我們的判斷條件有: 、
和
。我們分別定義:
綜合這三個條件就可以判斷在哪個扇區。那麼有沒有一種演算法可以將這一堆判斷數值化並轉換成1~6的數字呢?可以用下面的公式:
式中A代表 的正負,B代表 的正負,C代表 的正負,大於0為1,小於0為0。最後轉換出來的的N即為1~6的數字:
扇區
至此,我們成功完成了扇區判斷。
計算相鄰向量作用時長 控制相鄰向量作用時長就可以控制轉子到達任意方向,下面進行分析。
六個向量的大小 六個MOS管可以產生8種狀態
image-20220330152650462
設上開下合為0(電流從O往對應的向流),上合下開為1(電流從對應的嚮往O流),表示其中的六個向量。
image-20220330152634407
放在一張圖中即為:
image-20220330152741018
還有兩個零向量(000和111),無電流,不產生磁場。
對於100的狀態,可以等效為下面的電路圖:
image-20220330154344462
可以計算出馬達中三個相電壓(每相對於馬達中間連接點的電壓)
同理可以計算其他所有方向向量的相電壓,可以看出,六個向量的大小均為
,即SVPWM相電壓幅值為
電壓利用率 電壓利用率等於合成向量的電壓除以母線電壓。下面在複平面計算合成向量的電壓 :
根據歐拉公式可以推導出:
又因為三相電壓與相電壓幅值之間的關係:
帶入可以計算出 :
合成向量的電壓是相電壓幅值的
倍,而SVPWM相電壓幅值 為
,所以
即合成向量的電壓等於母線電壓。所以SVPWM的電壓利用率是100%。
SVPWM輸出電壓是馬鞍波 由於中間連接點N的點位是浮動的,為三角波,而相電壓是每相對於電機中間連接點N的電壓,所以相電壓不是一個正弦波,而是一個正弦波與一個三角波疊加而成的,即為馬鞍波。網圖:
d94c427b1d664ffe9b585d58e60482f6
向量作用時長 合成向量的電壓是所在磁區兩個向量與空向量不同時長的組合,其中 為空向量作用時長:
由於SVPWM的輸入是 和 ,但是要控制 和 ,所以要找出他們的對應關係。
對於第一個扇區,將 在 座標系中表示:
其中 和 根據前面的計算均為
,可以解出:
同理,可以計算出所有六個扇區:
第二區
第三區
第四扇區
第五扇區
第六扇區
六個扇區中都有相同的項,其中包含前文判斷扇區所用的 、 、 。直接把前面已經計算好的變數拿過來使用,大大減少了計算量。式中的
為調製比,定義:
可以將六個磁區表示為:
扇區 扇區 扇區 扇區 扇區 扇區
注意:當非零向量作用時間 ,需要進行過飽和處理:
定時器比較值計算 我們前文算出來的 和 以秒為單位,在單晶片中使用定時器控制MOS管的通斷需要配置比較值,所以需要把 和 轉換為三個互補的定時器比較值 、 和 。
定時器應設定為中心對齊模式,若為向上計數,計數器會從0開始計數到最大值,再反向從最大值計數到0;向下計數反之。故只需要控制前半個週期的比較值就可以產生相對中心對稱的PWM波。
CB0D861C81399F86BFB03E925F9AF8AC
為了減少MOS管的開關損耗,提高其使用壽命,應盡量減少MOS管的開關次數。在一個磁區內切換狀態的時候,合理使用零向量可以保證每次切換只改變一個MOS管。
先來看第一扇區,以每次只改變一個MOS管為原則,則切換順序為:
可以看到,切換順序構成了一個環路。在一個週期內我們需要控制三段作用時數:
六個MOS需要產生六路PWM來控制他們的狀態。由於上下半橋是互補的,所以只需要產生三個PWM:
image-20220331151522966
由於是中心對齊模式,所以只需要控制半個週期的時長
和
。在半個週期內 出現了兩次,分別為000和111,在半個週期內將這兩段時間平均分配。即為
。
這樣就可以計算出三個定時器的比較值 、 和 :
同理,可以計算出全部六個磁區的切換順序,這裡計算過程不再贅述。由影像和計算結果可以分析出,無論在哪個扇區中,三個定時器切換的圖像都是由上圖中三個方波構成,三個圖像排列組合也正好就是六個扇區的組合方式。故為了方便程式運行,將六個磁區的定時器比較值歸納如下:
設
此時六個扇區可以表示為:
扇區 扇區 扇區 扇區 扇區 扇區
也就是說程式設計時只需要計算 、 和 ,透過判斷就可以得到對應扇區的比較值。
至此,SVPWM輸出比較值 、 、 ,互補得到六個比較值,輸入到定時器,輸出三路PWM。
MATLAB仿真 假設需要產生10Khz的PWM波,則一個週期為0.0001秒。若單晶片主頻為180Mhz,預分頻係數為100-1,由下方計算公式:(其中 是計數值, 是預分頻值)
主頻
180,000,000/(18,000*100)=10,000Hz,則定時器的計數值應設定為18000-1。
在模型中設定 為24, 為18000。
image-20220331163919487
在foc子模組中,將park反變換的輸出輸入到SVPWM模組:
image-20220331163836635
unction [T1,T2,T3,sector] = fcn(Ualpha,Ubeta,Udc,Tpwm) % 初始化 sector = single(0); T1 = single(0); T2 = single(0); T3 = single(0); % 第一步:扇区判断 % 计算三个临时变量 Ua = Ubeta; Ub = (sqrt(3)*Ualpha - Ubeta)/2; Uc = (-sqrt(3)*Ualpha - Ubeta)/2; % 计算判断扇区所用的ABC的值 A=single(0);B=single(0);C=single(0);N=single(0); if(Ua>0) A = single(1); elseif(Ua<0) A = single(0); end if(Ub>0) B = single(1); elseif(Ub<0) B = single(0); end if(Uc>0) C = single(1); elseif(Uc<0) C = single(0); end % 计算判断扇区所用的N的值 N = A + 2*B + 4*C; % 扇区判断 switch (N) case 3 sector = single(1); case 1 sector = single(2); case 5 sector = single(3); case 4 sector = single(4); case 6 sector = single(5); case 2 sector = single(6); end % 第二步:计算相邻矢量作用时长 % 计算调制比 Umr = sqrt(3)*Tpwm/Udc; % 计算三个临时变量 X = Umr * Ua; Y = Umr * Ub; Z = Umr * Uc; % 分扇区计算Tx和Ty的值 Tx=single(0);Ty=single(0); switch (sector) case 1 Tx = Y; Ty = X; case 2 Tx = -Y; Ty = -Z; case 3 Tx = X; Ty = Z; case 4 Tx = -X; Ty = -Y; case 5 Tx = Z; Ty = Y; case 6 Tx = -Z; Ty = -X; end % 第三步:定时器比较值计算 % 干啥用的? if Tx+Ty > Tpwm Tx = Tx/(Tx+Ty); Ty = Ty/(Tx+Ty); else Tx = Tx; Ty = Ty; end % 计算三个临时变量 Ta = (Tpwm-Tx-Ty)/4.0; Tb = Ta+Tx/2.0; Tc = Tb+Ty/2.0; % 分扇区计算定时器比较值T1、T2和T3的值 switch (sector) case 1 T1 = Ta; T2 = Tb; T3 = Tc; case 2 T1 = Tb; T2 = Ta; T3 = Tc; case 3 T1 = Tc; T2 = Ta; T3 = Tb; case 4 T1 = Tc; T2 = Tb; T3 = Ta; case 5 T1 = Tb; T2 = Tc; T3 = Ta; case 6 T1 = Ta; T2 = Tc; T3 = Tb; end end
image-20220331165323898
扇區為從1到6的循環
image-20220331165348883
前文的狀態空間方程中 是非線性的,如果想對非線性系統進行卡爾曼濾波,需要對其線性化(Linearization),泰勒級數展開是將非線性系統線性化的方法。將 泰勒級數展開並取前兩項:
對 求偏導:
此時狀態空間方程式變為:
狀態方程式是建立在連續系統的基礎上,需要將其轉換成離散化系統後,才能在微控制器中實現卡爾曼濾波狀態觀測器。在離散化系統中,
,則狀態空間方程式表示為:
令 ,擴展卡爾曼濾波器的流程為以下五個步驟:
預測
計算預估值
計算誤差協方差
校正
計算卡爾曼增益
修正預估值
更新誤差協方差,用於下次計算
上面的五個方程中, 為採樣時間, 為模型誤差, 為測量誤差
function x_posteriori = fcn (Ialpha,Ibeta,Ualpha,Ubeta,Ls,Rs,Flux) global x; global P; R=[ 0.02 , 0 ; 0 , 0.02 ]; Q=[ 0.01 , 0 , 0 , 0 ; 0 , 0.01 , 0 , 0 ; 0 , 0 , 0.01 , 0 ; 0 , 0 , 0 , 0.010 , Time= 0 , Time ]; 0.0001 ; u=[Ualpha;Ubeta]; y=[Ialpha;Ibeta]; H=[ 1 , 0 , 0 , 0 ; 0 , 1 , 0 , 0 ]; B=[ 1 /Ls, 0 ; 0 , 1 /Ls ; 0 , 0 ; 0 , 0 ]; F=[-Rs/Ls, 0 ,Flux* sin (x( 4 , 1 ))/Ls,x( 3 , 1 )*Flux* cos (x( 4 , 1 ))/Ls; 0 ,-Rs/Ls,-Flux* cos (x( 4 , 1 ))/Ls,x( 3 , 1 )*Flux* sin (x( 4 , 1 ))/Ls; 0 , 0 , 0 , 0 ; 0 , 0 , 1 , 0 ]; f=[-Rs*Ialpha/Ls+x( 3 , 1 )*Flux* sin (x( 4 , 1 ))/Ls;-Rs*Ibeta /Ls-x( 3 , 1 )*Flux* cos (x( 4 , 1 ))/Ls; 0 ;x( 3 , 1 )]; % 計算預估值 x_prior=x+(f+B*u)* Time; % 計算誤差協方差 A=( eye ( 4 )+F*Time); P_prior=A*P*A'+Q; % 計算卡爾曼增益 K=(P_prior*H')/(H*P_prior* H'+R); % 修正預估值 x_posteriori=x_prior+K*(yH*x_prior); % 更新誤差協方差 P_posteriori=( eye ( 4 )-K*H)*P_prior; if x_posteriori( 4 , 1 )>( 2 * pi ) x_posteriori( 4 , 1 )=x_posteriori( 4 , 1 )-( 2 * pi ); end x=x_posteriori; P=P_posteriori;