s-function搭建变参数PMSM模型

S－function

一、什么是s-function

（1）非图形化地实现一个动态系统

（2）支持多种语言编写

（3）s-function是以现代控制理论的状态方程为基础的

二、如何编写一个自己的s-function

（1）介绍

支持m函数语言、C语言、C++语言、和Fortran

模板函数是s－function的框架，以下会介绍

S－function的参数有三个，基本只用到第一个“name”

（2）模板函数

只需要掌握（1）~（4）函数即可，（2）（3）分别用于连续和离散，一般都用离散的

t	时间
x	状态变量
u	输入
flag	s-function执行的标志，不用管
sys	输出
x0	状态变量初始值

第一部分：s-fucntion函数名为sfuntmpl

function [sys,x0,str,ts,simStateCompliance] = sfuntmpl(t,x,u,flag)
switch flag,
  case 0,   
[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;  	% （1）初始化函数，这个函数必须有

  case 1,
    sys=mdlDerivatives(t,x,u);                              % （2）连续下的微分输出，即dx
 
  case 2,
    sys=mdlUpdate(t,x,u);                                   % （3）离散下的微分输出，即x(k+1)
 
  case 3,
    sys=mdlOutputs(t,x,u);                                  % （4）输出sys
 
  case 4,
    sys=mdlGetTimeOfNextVarHit(t,x,u);						% （5）变 采样时间才会用到
 
  case 9,
    sys=mdlTerminate(t,x,u);
 
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
 
end

第二部分：初始化函数

初始化主要完成：

（1） 离散变量，输入输出个数的确定

（2） 状态变量的初始值

（3） 采样时间

function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes

sizes = simsizes;
 
sizes.NumContStates  = 0;           %连续状态变量的个数
sizes.NumDiscStates  = 0;           %离散状态变量的个数
sizes.NumOutputs     = 0;           %输出个数，
sizes.NumInputs      = 0;           %输入个数，

sizes.DirFeedthrough = 1;  			% 直馈通道，指输入在输出中直接显示出来

sizes.NumSampleTimes = 1;   		% at least one sample time is needed

sys = simsizes(sizes);

% initialize the initial conditions
%
x0  = [];               %状态变量的初始值
 
%
% str is always an empty matrix
%
str = [];
 
%
% initialize the array of sample times
%
ts  = [0 0];            %采样时间，ts = [Ts 0]

% Specify the block simStateCompliance. The allowed values are:
%    'UnknownSimState', < The default setting; warn and assume DefaultSimState
%    'DefaultSimState', < Same sim state as a built-in block
%    'HasNoSimState',   < No sim state
%    'DisallowSimState' < Error out when saving or restoring the model sim state
simStateCompliance = 'UnknownSimState';			% 不用管是什么

第三部分：连续微分函数

function sys=mdlDerivatives(t,x,u)
 
sys = [];

第四部分：离散差分函数

function sys=mdlUpdate(t,x,u)
 
sys = [];

第五部分：输出函数

function sys=mdlOutputs(t,x,u)
 
sys = [];

三、案例一：s-function来表示直流电压的动态模型

（1）Simulink中的直流电机的模型

这里为什么要加一个breaker，下面会讲到！

（2）用simulink搭建的直流电机模型

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Emg383XF-1621671565942)(file:///Z:/Users/Jia/AppData/Local/Temp/msohtmlclip1/01/clip_image002.jpg)]

电气元件都可以用电压/电流源的形式来表达

（3）s-function编写的直流电机模型

注意电流源的方向
$$U＝E_a+R{i}_a+L\frac{d{i}_a}{dt}$$
写成微分形式，
$$\frac{d{i}_a}{dt}=\left(u-K_e\omega -R{i}_a\right)/L$$

$$T_e–T_L=J\frac{d\omega }{dt}$$

$$J\frac{d\omega }{dt}=K_t{i}_a–T_L$$

根据状态方程，
$$\begin{gathered} \dot{x}=Ax+Bu \\ y=Cx+Du \end{gathered}$$

输入	端电压 u,负载转矩Tl
状态（可测量的,带微分的项）	转速w,和电枢电流ia
输出	也可以是转速w,和电枢电流ia

则，
$$\dot{x}=Ax+Bu$$

$$\dot{x}=\left[\begin{array}{cc}-\frac{R}{L}& -K\frac{e}{L}\\ K\frac{t}{J}& 0\end{array}\right]\left[\begin{array}{c}ia\\ w\end{array}\right]+\left[\begin{array}{cc}\frac{1}{L}& 0\\ 0& -\frac{1}{L}\end{array}\right]\left[\begin{array}{c}u\\ Tl\end{array}\right]$$

$$y=Cx$$

$$y=x=\left[\begin{array}{c}ia\\ w\end{array}\right]$$

用差分方程表示，
$$x\left(k+1\right)=x\left(k\right)+Ts\ast (Ax+Bu)$$
于是s-function为，

function [sys,x0,str,ts,simStateCompliance] = sfun_dc(t,x,u,flag)
 
switch flag,
 
  case 0,
    [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
 
  case 1,
    sys=mdlDerivatives(t,x,u);
 
  case 2,
    sys=mdlUpdate(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
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));
 
end
 
%====================================================================

function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
 
sizes = simsizes;
 
sizes.NumContStates  = 0;
sizes.NumDiscStates  = 2;
sizes.NumOutputs     = 2;
sizes.NumInputs      = 3;		%另外加了一个可变电阻，所以是三个输入变量
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;   % at least one sample time is needed
 
sys = simsizes(sizes);
 
%
% initialize the initial conditions
%
x0  = [0];
 
%
% str is always an empty matrix
%
str = [];
 
%
% initialize the array of sample times
%
Ts = 0.0001;
ts  = [Ts 0];
 
simStateCompliance = 'UnknownSimState';
 
%====================================================================

function sys=mdlDerivatives(t,x,u)
 
sys = [];
 
%====================================================================
function sys=mdlUpdate(t,x,u)
Ts=0.0001;				%这里要重新定义一次，变量不是全局的
 
 
Ke = 1.8*30/pi;
R = u(3);		%可变电阻
L = 0.02;
J = 10;
Kt = 17.2;
 
A = [-R/L -Ke/L; Kt/J 0];
B = [1/L 0; 0 -1/J];
 
sys = x+(A*x+B*[u(1);u(2)])*Ts;			% 此处的u要改成这个，不然维度不对
 
%====================================================================
function sys=mdlOutputs(t,x,u)
 
sys = x;
 
%====================================================================

function sys=mdlGetTimeOfNextVarHit(t,x,u)
 
sampleTime = 1;    %  Example, set the next hit to be one second later.
sys = t + sampleTime;
 
%====================================================================
function sys=mdlTerminate(t,x,u)
 
sys = [];
 
% end mdlTerminate

（4）仿真结果

原始模型不加breaker(在0s处合上)，不加的话，在0时刻以前就是电流值，所以下图的电流初始为5A。这是由于Laplace零时刻初始导致的。加上breaker后，变成了零状态，也就完全吻合了。

四、案例二：用s-function编写可变参数的PMSM模型

目的：

（1）电机参数不可在线动态修改

（2）simulink自带模型的坐标变换与常规逻辑不同

SIMULINK 中集成 PMSM 模块定义的坐标变换角与国内常规定义不同［4］，如图 1 所示。常规定义中 θ 为 d 轴与 α 轴的夹角，而 SIMULINK 中则定义成 q 轴与 α 轴的夹角，记为图中的 γ，显然，θ = γ －90°，这导致了在常规的调速系统中反馈回路需要对θ 做一次变换。

需求分析：

永磁同步电机的数学模型

（1）电气方程

显然，电气方程中的状态变量为id和iq，输入变量为ud和uq。但转速呢？在矩阵中，这算是输入还是状态？把它当作是输入量，但并不在u中体现。再增加一个输出为电磁转矩Te

function [sys,x0,str,ts,simStateCompliance] = Electrical_equation_signal_continus(t,x,u,flag)

switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,   
    [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;  %这个函数必须有

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u);                              %这个函数可以有

  %%%%%%%%%%
  % Update %
  %%%%%%%%%%
  case 2,
    sys=mdlUpdate(t,x,u);                                   %这个函数可以有

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u);                                  %这个函数必须有

  %%%%%%%%%%%%%%%%%%%%%%%
  % GetTimeOfNextVarHit %
  %%%%%%%%%%%%%%%%%%%%%%%
  case 4,
    sys=mdlGetTimeOfNextVarHit(t,x,u);

  %%%%%%%%%%%%%
  % Terminate %
  %%%%%%%%%%%%%
  case 9,
    sys=mdlTerminate(t,x,u);

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end

% end sfuntmpl

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes

%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded.  This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;

sizes.NumContStates  = 2;           %连续变量的个数，配置成连续
sizes.NumDiscStates  = 0;           %离散变量的个数
sizes.NumOutputs     = 2;           %输出个数，
sizes.NumInputs      = 8;           %输入个数，
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed

sys = simsizes(sizes);

%
% initialize the initial conditions
%
x0  = [0 0];               %状态变量的初始值

%
% str is always an empty matrix
%
str = [];

%
% initialize the array of sample times
%
ts  = [0 0];            



% Specify the block simStateCompliance. The allowed values are:
%    'UnknownSimState', < The default setting; warn and assume DefaultSimState
%    'DefaultSimState', < Same sim state as a built-in block
%    'HasNoSimState',   < No sim state
%    'DisallowSimState' < Error out when saving or restoring the model sim state
simStateCompliance = 'UnknownSimState';

% end mdlInitializeSizes

%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u)
                    % R = 2.875;
                    % Ld = 8.5e-3;
                    % Lq = 8.5e-3;
                    % Flux = 0.175;
we = u(3);
R = u(4);
Ld = u(5);
Lq = u(6);
Flux = u(7);

A = [   -R/Ld   we*Lq/Ld;
    -we*Ld/Lq   -R/Lq];
B = [    1/Ld    0;
            0 1/Lq];
C = [        0;
    -we*Flux/Lq];


sys = A*x + B*[u(1);u(2)] + C;


% end mdlDerivatives

%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
function sys=mdlUpdate(t,x,u)

sys = [];
% end mdlUpdate

%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u)

sys = [x];

% end mdlOutputs

%
%=============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block.  Note that the result is
% absolute time.  Note that this function is only used when you specify a
% variable discrete-time sample time [-2 0] in the sample time array in
% mdlInitializeSizes.
%=============================================================================
%
function sys=mdlGetTimeOfNextVarHit(t,x,u)

sampleTime = 1;    %  Example, set the next hit to be one second later.
sys = t + sampleTime;

% end mdlGetTimeOfNextVarHit

%
%=============================================================================
% mdlTerminate
% Perform any end of simulation tasks.
%=============================================================================
%
function sys=mdlTerminate(t,x,u)

sys = [];

% end mdlTerminate

（2）机械方程

function [sys,x0,str,ts,simStateCompliance] = Mechanical_equation_signal_continus(t,x,u,flag)

switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,   
    [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;  %这个函数必须有

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u);                              %这个函数可以有

  %%%%%%%%%%
  % Update %
  %%%%%%%%%%
  case 2,
    sys=mdlUpdate(t,x,u);                                   %这个函数可以有

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u);                                  %这个函数必须有

  %%%%%%%%%%%%%%%%%%%%%%%
  % GetTimeOfNextVarHit %
  %%%%%%%%%%%%%%%%%%%%%%%
  case 4,
    sys=mdlGetTimeOfNextVarHit(t,x,u);

  %%%%%%%%%%%%%
  % Terminate %
  %%%%%%%%%%%%%
  case 9,
    sys=mdlTerminate(t,x,u);

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end

% end sfuntmpl

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes

%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded.  This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;

sizes.NumContStates  = 1;           %连续变量的个数，配置成连续
sizes.NumDiscStates  = 0;           %离散变量的个数
sizes.NumOutputs     = 1;           %输出个数，
sizes.NumInputs      = 4;           %输入个数，
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed

sys = simsizes(sizes);

%
% initialize the initial conditions
%
x0  = [0];               %状态变量的初始值

%
% str is always an empty matrix
%
str = [];

%
% initialize the array of sample times
%
ts  = [0 0];            %采样时间，ts = [Ts 0]

% Specify the block simStateCompliance. The allowed values are:
%    'UnknownSimState', < The default setting; warn and assume DefaultSimState
%    'DefaultSimState', < Same sim state as a built-in block
%    'HasNoSimState',   < No sim state
%    'DisallowSimState' < Error out when saving or restoring the model sim state
simStateCompliance = 'UnknownSimState';

% end mdlInitializeSizes

%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u) 
Te = u(1);
TL = u(2);
J = u(3);
B = u(4);
% J = 0.001;
% B = 0.0008;

sys = (Te - TL - B*x)/J;


% end mdlDerivatives

%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
function sys=mdlUpdate(t,x,u)

sys = [];
% end mdlUpdate

%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u)

sys = x;

% end mdlOutputs

%
%=============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block.  Note that the result is
% absolute time.  Note that this function is only used when you specify a
% variable discrete-time sample time [-2 0] in the sample time array in
% mdlInitializeSizes.
%=============================================================================
%
function sys=mdlGetTimeOfNextVarHit(t,x,u)

sampleTime = 1;    %  Example, set the next hit to be one second later.
sys = t + sampleTime;

% end mdlGetTimeOfNextVarHit

%
%=============================================================================
% mdlTerminate
% Perform any end of simulation tasks.
%=============================================================================
%
function sys=mdlTerminate(t,x,u)

sys = [];

% end mdlTerminate

（3）用两个方程搭建信号模型

两个s-function分别见，

Electrical_equation_signal_continus.m

Mechanical_equation_signal_continus.m

注意以下几个点：

第一个坑：

做多参数输入的时候，电气方程不要输出转矩，会出现以下错误，维度没有错，但就是查不出原因。只能另外求

大概的原因是因为分母为零之类的。比如当我的J,B为0时，得到的错误信息是，

就是除了零导致的，但报的错误仍然是向量长度的问题。

网上的帖子：

输出Te后发现，果然是这个问题，计算结果出现NaN，一般有这样几种情况：0/0，Inf/Inf，Inf-Inf，Inf*0，这几种情况都会导致结果不确定，所以会得到NaN。

第二个坑：

因为反馈线id,iq搞反了，原本稳定的变乱了，查了一晚上，没结果，很影响心情。要及时备份，数据要检查，整理好。

（4）电路模型

在信号模型的基础上，再封装一层即可。

第一个坑：

由于自定义的 PMSM 模型属于信号模型，省去了功率部分，故可由正弦信号直接驱动，这就意味着连SVPWM和逆变器都不需要了。但是这好像不合理呀。配合适当的信号增益可表示功率大小。当控制信号采用 SPWM算法时，调制函数为正弦波，可以直接驱动电机模型。对于 SVPWM 算法，相电压的调制信号是马鞍形波［5］，不可直接引入电机模型，由于其线电压调制函数是正弦波，故可利用线电压信号驱动自定义模型. ，此时坐标变换[1]需改写成

---

[1] 参考文献《SIMULINK中PMSM模型的改进及在参数辨识中的应用_王莉娜》

第二个坑：

我们所搭建的pmsm模型是其数学模型，并不是电路模型，其端口电压是通过测量后再输入的。如何人为地将其进行再加工，转换成电路模型呢？

答：用受控电流源将信号模型变为电路模型

电路模型验证：

仿真文件见，Result04_closeloop_inverter_svpwm.slx

声明：案例一使用了B站UP主陈诚电气部分，如有侵权，请及时联系
仿真文件获取：https://download.csdn.net/download/weixin_42362056/18970461