如何解决simulink控制系统仿真中的代数环

1. 什么是代数环

2. 如何解决代数环

3. 多个s函数导致的代数环

4. 源代码

1. 什么是代数环

在simulink仿真过程中，当输入信号直接取决于输出信号，同时输出信号也直接取决于输入信号时，由于数字计算的时序性，而出现的由于没有输入无法计算输出，没有输出也无法得到输入的“死循环” ，称之为代数环。

如下图所示，output = func（input+output）。初始时，由于没有output，所以不能计算func函数；但是为了得到output，又必须要计算func。如此往复就形成了代数环。

2. 如何解决代数环

——连续模型

对于连续模型，可以在反馈通道添加Memory模块延时

——离散模型

对于离散模型，可以在反馈通道添加delay模块延时

3. 多个s函数导致的代数环

对于控制系统的仿真来说，很少会出现代数环的问题。下图这种系统就会出现代数环，但是这种系统本身就不存在（至少控制系统不会出现这种玩意）。

当simulink中的模块不足以完成建模时，常用的方法就是使用s函数了。如果只是一个s函数建立的模型，例如下图，也不会出现代数环的问题（如果简单粗暴的用s函数建立上图那种模型当我没说）。

如果用到多个s函数，如下图，就会提示代数环的警告了。

systemP的模型：（源代码放最后啦）

$\begin{bmatrix} \dot{x}1(t)\\ \dot{x}2(t) \end{bmatrix} = \begin{bmatrix} 0 & 1\\ 0 & -1 \end{bmatrix}\begin{bmatrix} x1(t))\\ x2(t)) \end{bmatrix}+\begin{bmatrix} 0\\1 \end{bmatrix} u(t) - \begin{bmatrix} 0\\-\frac{1}{2} \end{bmatrix} M(t)$

$y = \begin{bmatrix} x1(t)\\ x2(t)\end{bmatrix}$

systemM的模型：（源代码放最后啦）

$M1 = 19.98 * sign(n0 - n) * (0.005 * sin(a) + 0.02 * cos(a))$ ①

$Mmud = 2 / 15 * pi^2 * 0.05 * 1 * (0.0675^2 * 0.055^2 / (0.0675^2 - 0.055^2)) * (n0 - n)$ ②

$Mw=1.7987 * (tanh(32.8317 * w)-tanh(7.4921 * w)) + 0.9601 * tanh(6.0375 * w) + 0.0378 * w$ ③

$Mf = 0.1 * sin(w) + 0.32$ ④

$y = M1+Mmud+Mw+Mf$

对于这种情况如果我们加入Memory模块做延时，虽然能解决代数环的问题，但是如果系统本身不是稳定的，延时环节会让系统产生震荡。

将加入Memory和没加Memory的两种输出的波形图对比：（本身搭建的系统就没有意义，所以不分析，将分析另外一个波形图，以此说明差异）

上面的两个模型没有严密的逻辑分析，紧供分析代数环，所以对于系统震荡的问题不好分析，所以在此插入我另外一个仿真中遇到的两个图（也是systemM和systemP构成的仿真，只是系统严谨）

下图是加入Memory延时，导致机器学习训练收敛性很差

下图是去掉延时，解决代数环后，通过机器学习后的控制波形

很容易发现，第二个的效果很好啦

延时对不稳定系统的影响还是很大的。对于多个s函数，导致计算时序的影响，出现代数环的问题，我们可以在s函数上做调整。

在s函数，flag等于0时，也就是初始化时，我们可以看到

sizes.DirFeedthrough

的设置，输入是否直接影响输出。将其值设为0(修改连续系统的DirFeedthrough，如果修改纯数学公式的会导致仿真错误，本文中systemP为连续系统，systemM为数学公式)；仿真时就不会出现代数环的问题，且波形和有代数环时一样。

我在做仿真的时候，systemP和systemM可以用simulink提供的模块搭建，但是会比较麻烦，所以采用了s函数，后来也做了完全用simulink提供的模块搭建的仿真，控制效果也和s函数出现的效果类似。

4. 源代码

systemP：

%systemP
function [sys,x0,str,ts,simStateCompliance] = systemP(t,x,u,flag)switch flagcase 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));
endfunction [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates  = 2;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 2;
sizes.NumInputs      = 2;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;  sys = simsizes(sizes);x0  = [0;0];str = [];ts  = [0 0]; simStateCompliance = 'UnknownSimState';function sys=mdlDerivatives(t,x,u)
sys = [0 1;0 -1] * x + [0;1]*u(1) - [0;-1/2]*u(2);function sys=mdlUpdate(t,x,u)
sys = [];function sys=mdlOutputs(t,x,u)
sys = x;function sys=mdlGetTimeOfNextVarHit(t,x,u)
sampleTime = 1;
sys = t + sampleTime;function sys=mdlTerminate(t,x,u)
sys = [];

systemM：

%systemM
function [sys,x0,str,ts,simStateCompliance] = systemM(t,x,u,flag)switch flagcase 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));
endfunction [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates  = 0;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 1;
sizes.NumInputs      = 3;
sizes.DirFeedthrough = 1;
sizes.NumSampleTimes = 1;  sys = simsizes(sizes);x0  = [];str = [];ts  = [0 0]; %[0 0] simStateCompliance = 'UnknownSimState'; function sys=mdlDerivatives(t,x,u)
sys = [];function sys=mdlUpdate(t,x,u)
sys = [];function sys=mdlOutputs(t,x,u)
n0 = u(1);        %r/min
a = u(2)*pi/180;  %degree
w = u(3);         %rad/s
n = w*60/(2*pi);  %r/minM1 = 19.98 * sign(n0 - n) * (0.005 * sin(a) + 0.02 * cos(a));
Mmud = 2 / 15 * pi^2 * 0.05 * 1 * (0.0675^2 * 0.055^2 / (0.0675^2 - 0.055^2)) * (n0 - n);
Mw = 1.7987 * (tanh(32.8317 * w)-tanh(7.4921 * w)) + 0.9601 * tanh(6.0375 * w) + 0.0378 * w;
Mf = 0.1 * sin(w) + 0.32;
ML = M1+Mmud+Mw+Mf;
sys = ML;function sys=mdlGetTimeOfNextVarHit(t,x,u)
sampleTime = 1;
sys = t + sampleTime;function sys=mdlTerminate(t,x,u)
sys = [];