matlab中S函数的概念及使用
??S函數即系統函數System Function的意思，為什么要使用S函數呢？是因為在研究中，有時需要用到復雜的算法設計等，而這些算法因為其復雜性不適合用普通的Simulink模塊來搭建，即matlab所提供的Simulink模塊不能滿足用戶的需求，需要用編程的形式設計出S函數模塊，將其嵌入到系統中。如果恰當地使用S函數，理論上，可以在Simulink下對任意復雜的系統進行仿真。
??? S函數具有固定的程序格式，用matlab語言可以編寫S函數，此外還允許用戶使用C、C++、Fortran和Ada等語言進行編寫，用非matlab語言進行編寫時，需要采用編譯器生成動態鏈接庫DLL文件。
在主窗口中輸入sfundemos，或者點擊Simulink->User-Defined Functions->S-Function Examples，即可出現如圖1所示的界面，可以選擇對應的編程語言查看演示文件。
圖1 S函數范例庫
Matlab為了用戶使用方便，有一個S函數的模板sfuntmpl.m，一般來說，我們僅需要在sfuntmpl.m的基礎上進行修改即可。在主窗口輸入edit sfuntmpl即可出現模板函數的內容，可以詳細地觀察其幫助說明以便更好地了解S函數的工作原理。模板函數的定義形式為function[sys,x0,str,ts]=sfuntmpl(t,x,u,flag)，一般來說，S函數的定義形式為[sys,x0,str,ts]=sfunc(t,x,u,flag,p1,…Pn)，其中的sfunc為自己定義的函數名稱，以上參數中，t、x、u分別對應時間、狀態、輸入信號，flag為標志位，其取值不同，S函數執行的任務和返回的數據也是不同的，pn為額外的參數，sys為一個通用的返回參數值，其數值根據flag的不同而不同，x0為狀態初始數值，str在目前為止的matlab版本中并沒有什么作用，一般str=[]即可，ts為一個兩列的矩陣，包含采樣時間和偏移量兩個參數，如果設置為[0 0]，那么每個連續的采樣時間步都運行，[-1 0]則表示按照所連接的模塊的采樣速率進行，[0.25 0.1]表示仿真開始的0.1s后每0.25s運行一次，采樣時間點為TimeHit=n*period+offset。
S函數的使用過程中有2個概念值得注意：1、direct feedthrough，系統的輸出是否直接和輸入相關聯，即輸入是否出現在輸出端的標志，若是為1，否則為0，一般可以根據在flag=3的時候，mdlOutputs函數是否調用輸入u來判斷是否直接饋通。2、dynamically sized inputs，主要給出連續狀態的個數、離散狀態的個數、輸入數目、輸出數目和直接饋通否。
S函數中目前支持的flag選擇有0、1、2、3、4、9等幾個數值，下面說一下在不同的flag情況下S函數的執行情況。1）flag=0。進行系統的初始化過程，調用mdlInitializeSizes函數，對參數進行初始化設置，比如離散狀態個數、連續狀態個數、模塊輸入和輸出的路數、模塊的采樣周期個數、狀態變量初始數值等。2）flag=1。進行連續狀態變量的更新，調用mdlDerivatives函數。3）flag=2。進行離散狀態變量的更新，調用mdlUpdate函數。4）flag=3。求取系統的輸出信號，調用mdlOutputs函數。5）flag=4。調用mdlGetTimeOfNextVarHit函數，計算下一仿真時刻，由sys返回。6）flag=9。終止仿真過程，調用mdlTerminate函數。
圖2 不同flag情況下S函數執行情況
在實際仿真過程中，Simulink會自動將flag設置為0，進行初始化過程，然后將flag的數值設置為3，計算模塊的輸出，一個仿真周期后，Simulink將flag的數值先后設置為1和2，更新系統的連續和離散狀態，再將其設置為3，計算模塊的輸出，如此循環直至仿真結束條件滿足。?
在S函數的編寫過程中，首先需要搞清楚模塊中有多少個連續和離散狀態，離散模塊的采樣周期是如何的，同時需要了解模塊的連續和離散的狀態方程分別是什么，輸出如何表示。下面以實例說明S函數的具體應用。
實例一：在進行BLDC的速度電流雙閉環控制時候，如果電流采用滯緩控制，那么速度環的輸出作為參考電流的輸入Is，A、B、C三相的參考電流分別為I_ar、I_br、I_cr，轉子的位置記為Pos，那么轉子位置和三相參考電流之間的關系表如圖3所示。
圖3 轉子位置和三相參考電流的關系表
分析系統情況：輸入為兩路，一是角度Angle，通過mod(angle,2*pi)轉化為Pos，二是參考電流輸入Is，輸出為三路，分別為A、B、C的三相參考電流，離散和連續的狀態均為0，輸入端出現在輸出端，將S函數命名為given_current，其程序代碼如下所示：
%參考電流模塊
%Author:dingqian12345@126.com
%很明顯：2個輸入分別為角度和速度控制器的輸出、3個輸出為三相電流的參考電流、為直通模型
function [sys,x0,str,ts]=cemf(t,x,u,flag)
clc;
switch flag
??? case 0
??????? [sys,x0,str,ts]=mdlInitializeSizes;???? %初始化
??? case 1
??????? sys=[];???????????????????????????????? %連續狀態的更新
??? case 2
??????? sys=[];????????????? %離散狀態的更新
??? case 3
??????? sys=mdlOutputs(u);??????????? ??????????%求取系統的輸出信號
??? case 4
??????? sys=[];???????????????????????????????? %計算下一時刻的仿真時間
??? case 9
??????? sys=[];???????????????????????????????? %終止仿真
??? otherwise
??????? error(['Unhandled flag=',num2str(flag)]);
end
%%%%%在flag=0的時候進行整個系統的初始化
function [sys,x0,str,ts]=mdlInitializeSizes
sizes = simsizes;??????????? %讀入初始化參數模板
sizes.NumContStates? = 0;??? %連續狀態個數
sizes.NumDiscStates? = 0;??? %離散狀態個數
sizes.NumOutputs???? = 3;??? %輸出變量個數
sizes.NumInputs????? = 2;??? %輸入信號個數
sizes.DirFeedthrough = 1;? ??%輸入直接傳入輸出信號否
sizes.NumSampleTimes = 1;??? % at least one sample time is neededWO? 一般來說為1個
sys = simsizes(sizes);
x0=[];??? %狀態初始化
str=[];
ts=[-1 0];??? %采樣周期若寫成-1表示繼承其輸入信號
%%%%%%在flag=1的時候進行連續系統狀態的更新
function sys=mdlDerivatives(t,x,u)
sys = [];
%%%%%%在flag=2的時候進行離散系統狀態的更新
function sys=mdlUpdates
%sys(1,1)=x(1)+T*x(2);??????????????? %為什么會寫成這樣呢
%sys(2,1)=x(2)+T*fst2(x,u,r,h);
%%%%%%在flag=3的時候進行系統輸出信號的求取
function sys=mdlOutputs(u)
m=current_dq(u(1),u(2));
sys(1,1)=m(1);
sys(2,1)=m(2);
sys(3,1)=m(3);
%%%%%%在flag=4的時候進行下一時刻仿真時間的計算
function sys=mdlGetTimeOfNextVarHit(t,x,u)
%sampleTime = 1;??? %? Example, set the next hit to be one second later.
%sys = t + sampleTime;
%%%%%%在flag=9的時候終止仿真過程
function sys=mdlTerminate(t,x,u)
%sys = [];
??? 如果某個flag對應的函數不起作用，我們將其所調用的函數設為空或者sys=[]，為了保持模塊化以便思路清晰和方便擴充，一般不要使用case {1,4,9 }的形式，而是一個一個的列出。在以上代碼中，在flag=3的時候調用了自編的子程序m=current_dq(u(1),u(2))，current_dq的代碼如下：
%%%%%%用戶自定義子程序
function x=current_dq(angle,current)
pos=mod(angle,2*pi);
x(1)=current;
x(2)=-current;
x(3)=0;
if 0<=pos && pos<pi/3
???? x(1)=current;
???? x(2)=-current;
???? x(3)=0;
end
if pi/3<=pos && pos<2*pi/3
???? x(1)=current;
???? x(2)=0;
???? x(3)=-current;
end
if 2*pi/3<=pos && pos<pi
??? x(1)=0;
??? x(2)=current;
??? x(3)=-current;
end
if pi<=pos && pos<4*pi/3
??? x(1)=-current;
??? x(2)=current;
??? x(3)=0;
end
if 4*pi/3<=pos && pos<5*pi/3
??? x(1)=-current;
??? x(2)=0;
??? x(3)=current;
end
if 5*pi/3<=pos && pos<2*pi
?? x(1)=0;
?? x(2)=-current;
?? x(3)=current;
end
??? S函數的代碼完成后，下面進行Simulink模塊的搭建，如圖4所示，其中的S-Function模塊提取路徑為Simulink->User-Defined Functions->S-Function，雙擊S-Function模塊，在name行中輸入函數名稱given_current即可，仿真結果如圖5所示。
圖4 參考電流模塊Simulink測試圖
圖5 參考電流模塊測試結果
????實例二：我覺得對于不做理論研究的我來說，如果用到S函數的話，一般的情況就如實例一的情況那樣，沒有離散和連續的狀態，只需要對輸入做某些運算或者判斷就可。上學期學的一門交流調速的課程，問題稍微復雜一點，在此也說一下。
??? 對于含有齒隙環節的閉環系統，系統控制量反映了系統運動方向的變化，可以對系統控制量U加一個補償量Ucomp以縮短系統在齒隙內的運行時間，補償控制系統框圖如圖6所示。具體的補償規則為：
if U(k)>0 and U(k-1)<=0 then flag=1
if U(k)<0 and U(k-1)>=0 then flag=1 Ucomp=flag*const*exp(-t/T0)
圖6 補償控制系統框圖
??? 分析系統情況：需離散狀態變量，系統狀態方程X1(k+1)=U(k)、X2(k+1)=f(X1(k)，U(k))，其中f代表一個如上等式所示的關系，輸出Y= flag*const*exp(-t/T0)，因此得到如下的信息，輸入為一路，即為圖6中的u，輸出為一路，連續狀態個數為0，離散狀態個數為2，輸入端不出現在輸出端，將S函數命名為comp，根據以上分析，其程序代碼如下所示：
%Author:dingqian12345@126.com
%Note:實現基于規則的補償算法
%2個離散狀態、0個連續狀態、1個輸入信號、1個輸出信號、輸出不涉及輸入信號
function [sys,x0,str,ts]=comp(t,x,u,flag)
clc;
global Cflag;
global Cont;
global T0;
Cont=20;
T0=1.5;
switch flag
??? case 0
??????? [sys,x0,str,ts]=mdlInitializeSizes; %初始化
??? case 1
??????? sys=mdlDerivatives(t,x,u);????????? %連續狀態的更新
??? case 2
??????? sys=mdlUpdates(t,x,u);???? ?????????%離散狀態的更新
??? case 3
??????? sys=mdlOutputs(t,x,u);????????????????? %求取系統的輸出信號
??? case 4
??????? sys=mdlGetTimeOfNextVarHit(t,x,u);????? %計算下一時刻的仿真時間
??? case 9
??????? sys=mdlTerminate(t,x,u);??????????? %終止仿真
??? otherwise
??????? error(['Unhandled flag=',num2str(flag)]);
end
%%%%%在flag=0的時候進行整個系統的初始化
function [sys,x0,str,ts]=mdlInitializeSizes
sizes = simsizes;??????????? %讀入初始化參數模板
sizes.NumContStates? = 0;??? %連續狀態個數
sizes.NumDiscStates? = 2;??? %離散狀態個數
sizes.NumOutputs???? = 1;??? %輸出變量個數
sizes.NumInputs????? = 1;??? %輸入信號個數
sizes.DirFeedthrough = 0;??? %輸入直接傳入輸出信號否
sizes.NumSampleTimes = 1;??? % at least one sample time is neededWO? 一般來說為1個
sys = simsizes(sizes);
x0=[0;0];??????? %狀態初始化
str=[];
ts=[-1 0];??? %采樣周期若寫成-1表示繼承其輸入信號
%%%%%%在flag=1的時候進行連續系統狀態的更新
function sys=mdlDerivatives(t,x,u)
sys = [];
%%%%%%在flag=2的時候進行離散系統狀態的更新
function sys=mdlUpdates(t,x,u)
global Cflag;
global Cont;
global T0;
Cflag=comp_decide(t,x,u);
sys(1,1)=u;
sys(2,1)=Cflag*Cont*exp(-t/T0);
%%%%%%在flag=3的時候進行系統輸出信號的求取
function sys=mdlOutputs(t,x,u)
sys=x(2);
%%%%%%在flag=4的時候進行下一時刻仿真時間的計算
function sys=mdlGetTimeOfNextVarHit(t,x,u)
%sampleTime = 1;??? %? Example, set the next hit to be one second later.
%sys = t + sampleTime;
%%%%%%在flag=9的時候終止仿真過程
function sys=mdlTerminate(t,x,u)
sys = [];
??? 可以看出，這里的S函數比實例1的S函數多用到了離散狀態更新函數mdlUpdates，注意狀態的更新方程的寫法，sys=x(n+1)，可以通過help sfuntmpl查看。在以上代碼中，在flag=2的時候調用了自編的子程序m=comp_decide(t,x,u)，comp_decide的代碼如下：
%%%%%%用戶自定義子程序
function cc=comp_decide(t,x,u)
cc=0;
if u>0 && x(1)<=0
??? cc=1;
end;
if u<0 && x(1)>=0
??? cc=-1;
end;
??? 最后，將S函數嵌入到Simulink中，所搭建的補償仿真框圖如圖7所示，我們通過Scope4觀察控制量U和補償Ucomp，其中cont為20，To=1.5。圖8為補償系統中U和Ucomp仿真結果局部放大圖，可見，系統可以捕捉到控制量過零點的時刻并給出相應的補償量。
圖7 補償系統Simulink仿真框圖
圖8 補償系統中U和Ucomp仿真結果局部放大圖
?
參考文獻
[1] 薛定宇 著.控制系統仿真與計算機輔助設計[M].北京：機械工業出版社出版社，2009.
[2] 很好的S函數學習資料-Simulink的S-Function編寫指導[EB/OL].
http://www.ilovematlab.cn/thread-71226-1-1.html.[2010-12-26].
[3] Simulink-s_function使用及應用舉例[EB/OL].
http://www.ilovematlab.cn/thread-62299-1-1.html.[2010-12-26].
[4] 祁麟.交流伺服系統齒隙非線性控制及其網絡化研究[D].碩士論文.南京理工大學，2006.
?
2010年12月26日晚上完成于njust 電工樓 310房間
?
CopyRight：版權所有若需轉載或使用請聯系作者
Email:dingqian12345@126.com
總結
以上是生活随笔為你收集整理的matlab中S函数的概念及使用的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： mysql 设置事物自动提交_mysql
- 下一篇： matlab自带kfcm函数,kfcmF