目录

 1. 什么是代数环

  2. 如何解决代数环

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

  4. 源代码

1. 什么是代数环

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

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_13,color_FFFFFF,t_70,g_se,x_16

2. 如何解决代数环

——连续模型

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

5b032f0514374f0581dceda54ce3474d.png                   watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

——离散模型

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

d9129949b3164807afedfe68c2ff86e2.png                watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

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

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_11,color_FFFFFF,t_70,g_se,x_16

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_17,color_FFFFFF,t_70,g_se,x_16

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

systemP的模型:(源代码放最后啦)

gif.latex?%5Cbegin%7Bbmatrix%7D%20%5Cdot%7Bx%7D1%28t%29%5C%5C%20%5Cdot%7Bx%7D2%28t%29%20%5Cend%7Bbmatrix%7D%20%3D%20%5Cbegin%7Bbmatrix%7D%200%20%26%201%5C%5C%200%20%26%20-1%20%5Cend%7Bbmatrix%7D%5Cbegin%7Bbmatrix%7D%20x1%28t%29%29%5C%5C%20x2%28t%29%29%20%5Cend%7Bbmatrix%7D+%5Cbegin%7Bbmatrix%7D%200%5C%5C1%20%5Cend%7Bbmatrix%7D%20u%28t%29%20-%20%5Cbegin%7Bbmatrix%7D%200%5C%5C-%5Cfrac%7B1%7D%7B2%7D%20%5Cend%7Bbmatrix%7D%20M%28t%29

gif.latex?y%20%3D%20%5Cbegin%7Bbmatrix%7D%20x1%28t%29%5C%5C%20x2%28t%29%5Cend%7Bbmatrix%7D

systemM的模型:(源代码放最后啦)

gif.latex?M1%20%3D%2019.98%20*%20sign%28n0%20-%20n%29%20*%20%280.005%20*%20sin%28a%29%20+%200.02%20*%20cos%28a%29%29       ①

%20%280.0675%5E2%20-%200.055%5E2%29%29%20*%20%28n0%20-%20n%29   ②

 gif.latex?Mw%3D1.7987%20*%20%28tanh%2832.8317%20*%20w%29-tanh%287.4921%20*%20w%29%29%20+%200.9601%20*%20tanh%286.0375%20*%20w%29%20+%200.0378%20*%20w    ③

gif.latex?Mf%20%3D%200.1%20*%20sin%28w%29%20+%200.32       ④

gif.latex?y%20%3D%20M1+Mmud+Mw+Mf

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16


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

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

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


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

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

sizes.DirFeedthrough

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

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16

watermark,type_d3F5LXplbmhlaQ,shadow_50,text_Q1NETiBAU3Vkb1JlYm9vdA==,size_20,color_FFFFFF,t_70,g_se,x_16


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


4. 源代码

systemP:

%systemP
function [sys,x0,str,ts,simStateCompliance] = systemP(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  = 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 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  = 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/min

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; 
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 = [];

Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐