2020-02-23-Cassie电弧模型离散化

Posted on | In

对于Cassie电弧方程的离散化。

Cassie电弧方程为

\[\frac{1}{g} \cdot \frac{\mathrm{d} g}{\mathrm{d} t}=\frac{1}{\tau_{\mathrm{c}}}\left(\frac{e^{2}}{e_{0}^{2}}-1\right)\]

可以将$1/g$进行积分后，写入等式为

\[\frac{ { {\text{d} }\ln g} } { {dt} } = \frac{1} {\tau }(\frac{ { {u^2} } }{ { U_c^2} } - 1)\]

然后对右侧积分求出$lng$的值。

在Output子方法中，用$g=e^{lng}$来求得。

那么离散化只需要对右侧的积分进行离散化即可。

如果不提高仿真精度，那么等式右侧就是一个 常数，那么它的积分就是$lng = A * t$

那么用模型对比看一看效果吧。

这样做会报错，还需要进一步检查问题在哪。

上面这种离散化的方式不对。

按照微分写成离散化可以写成：

$\ln g ( t + 1 ) = \frac { 1 } { \tau_c } \left( \frac { e ^ { 2 }( t ) } { e ^ { 2 }_0 } - 1 \right) \cdot \Delta t + \ln g (t)$

那么，Update子函数应该怎么用呢？

在Update子函数中，是否可以利用Dwork来设置两个个全局变量，然后每次进行更新？

需要试验一下，Dwork的继承效果，先建一个简单的程序试一下，使之输出一个方波。

程序如下

function breaker(block)
% smileland作品，免费使用，但请不要删除此句话
% 以下设定参数并不一定全部需要，不需要的注释掉就好
setup(block);

%endfunction

function setup(block)
  %block.CurrentTime 只读对象，调用当前时间

  % 初始化对话框参数，即对话框内给定的初始化参数.
  block.NumDialogPrms     = 1;%给定3个参数
  %block.DialogPrmsTunable = {'Tunable','Nontunable','SimOnlyTunable'};%参数是否可变（此句话可选是否需要，默认不可变）

  %%%%%%%%初始化输入输出端口的属性%%%%%%%
  
  % 初始化输入输出端口个数.
  block.NumInputPorts  = 1;%输入端口个数
  block.NumOutputPorts = 1;%输出端口个数
  
  % 设定输入输出端口时继承还是动态，以指示输入和输出端口从模型继承其已编译属性（尺寸、数据类型、复杂度和采样模式）.
  block.SetPreCompInpPortInfoToDynamic;%设定输入端口为动态
  block.SetPreCompOutPortInfoToDynamic;%设定输出端口为动态

  % 设定输入端口属性.（也可以再Dwork中设定，在PostPropagationSetup子方法中）
  block.InputPort(1).DatatypeID  = 0;  % （-1 inherited；0 double;1 single;2 int8;3 uint8;4 int16;5 uint16;6 int32;7 uint12;8 boolean）
  block.InputPort(1).Complexity  = 'Real';%是否为复数
 
  % 设定输出端口属性.（也可以再Dwork中设定，在PostPropagationSetup子方法中）
  block.OutputPort(1).DatatypeID  = 0; % double
  block.OutputPort(1).Complexity  = 'Real';
  
  block.InputPort(1).Dimensions = 1;%输入口1的维数是1，最高支持到2维
  % block.InputPort(1).DimensionsMode = 'Fixed';%Fixed表示端口维数固定，Variable表示端口维数可变
  block.OutputPort(1).Dimensions = 1;%输出口1的维数是1，最高支持到2维
  % block.OutputPort(1).DimensionsMode = 'Fixed';%Fixed表示端口维数固定，Variable表示端口维数可变
  
  
  block.InputPort(1).DirectFeedthrough = true;%输入口有无直接馈通（true/false）

  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % 设定有多少连续状态变量.
  %block.NumContStates = 1;
  
  % 设定有多少离散状态变量.
  % block.NumDworkDiscStates = 1;

  % 设定采样时间.
  %  [0 0.5]            : 从0.5s开始连续采样
  %  [0.1 0.5]          : 从0.5s开始，每0.1s采样一次
  %  [-1, 0]            : 继承采样频率
  %  [-2, 0]            : 动态变化的采样频率
  block.SampleTimes = [2.5e-5 0];%连续采样率
  
  %设定离散工作向量个数
  % block.NumDworks = 1;
  
  % -----------------------------------------------------------------
  % Options，若使用TLC则需要以下语句，否则注释掉即可
  % -----------------------------------------------------------------

  % block.SetAccelRunOnTLC(false);
  
  % block.OperatingPointCompliance = 'Default';
  
  % -----------------------------------------------------------------
  % 引用使用的子方法，主要有PostPropagationSetup;InitializeConditions;
  % Start;Outputs;Updates;Derivatives;Terminate.以下子方法仅保留着几个。
  % -----------------------------------------------------------------
  
  %
  %包含PostPropagationSetup子方法，用于初始化DWork工作向量
  block.RegBlockMethod('PostPropagationSetup', @DoPostPropSetup);
  %}
  
  %包含InitializeConditions子方法，用于初始化状态变量或DWork工作向量的值
  block.RegBlockMethod('InitializeConditions', @InitializeConditions);
  
  %与InitializeConditions子方法功能一致，但只初始化一次，Ini方法每次调用子方法时均会初始化
  % block.RegBlockMethod('Start', @Start);

  %包含输出模块
  block.RegBlockMethod('Outputs', @Outputs);

  %包含离散状态更新模块
  block.RegBlockMethod('Update', @Update);

  %包含连续状态微分模块
  %block.RegBlockMethod('Derivatives', @Derivatives);
  
  %包含终止模块
  %block.RegBlockMethod('Terminate', @Terminate);

  %还有其他模块，仅将最主要的列出，若需要则使用S函数自带模板。
%endfunction

% -------------------------------------------------------------------
% 以下为子方法模块，同样主要有DoPostPropSetup;InitializeConditions;
  % Start;Outputs;Updates;Derivatives;Terminate.以下子方法仅保留这几个。
% -------------------------------------------------------------------
%
function DoPostPropSetup(block)
  block.NumDworks = 4;
  
  block.Dwork(1).Name            = 'out';
  block.Dwork(1).Dimensions      = 1;
  block.Dwork(1).DatatypeID      = 0;      % double
  block.Dwork(1).Complexity      = 'Real'; % real
  block.Dwork(1).UsedAsDiscState = true;
  
  block.Dwork(2).Name            = 'out0';
  block.Dwork(2).Dimensions      = 1;
  block.Dwork(2).DatatypeID      = 0;      % double
  block.Dwork(2).Complexity      = 'Real'; % real
  block.Dwork(2).UsedAsDiscState = true;
    
  block.Dwork(3).Name            = 'LastTime';
  block.Dwork(3).Dimensions      = 1;
  block.Dwork(3).DatatypeID      = 0;      % double
  block.Dwork(3).Complexity      = 'Real'; % real
  block.Dwork(3).UsedAsDiscState = true;
  
  block.Dwork(4).Name            = 'stage';
  block.Dwork(4).Dimensions      = 1;
  block.Dwork(4).DatatypeID      = 0;      % double
  block.Dwork(4).Complexity      = 'Real'; % real
  block.Dwork(4).UsedAsDiscState = true;
  
  % Register all tunable parameters as runtime parameters.
  block.AutoRegRuntimePrms;

%endfunction
%}
%

function InitializeConditions(block)

block.Dwork(1).Data = 0;
block.Dwork(3).Data = 0;

    
%{
function Start(block)

  tau = block.DialogPrm(1).Data;
  Uc = block.DialogPrm(2).Data; 

%}

%
  function Update(block)

      delta = block.DialogPrm(1).Data;
      
      if (block.CurrentTime - block.Dwork(3).Data) <= delta
          block.Dwork(1).Data = 0;
      elseif ((block.CurrentTime - block.Dwork(3).Data) > delta) && ((block.CurrentTime - block.Dwork(3).Data) < 2*delta)
          block.Dwork(1).Data = 1;
      else
          block.Dwork(3).Data = block.CurrentTime;
      end
      
%}

function Outputs(block)
  
  block.OutputPort(1).Data = block.InputPort(1).data * block.Dwork(1).Data;
  
%{
function Derivatives(block)

tau = block.DialogPrm(1).Data;
Uc = block.DialogPrm(2).Data; 

if block.CurrentTime <= block.DialogPrm(3).Data
    block.Dwork(3).Data = 0;
else
    block.Dwork(3).Data = 1;
end
    
block.Derivatives.Data = block.Dwork(3).Data * 1/tau * (block.InputPort(1).data^2 / Uc^2 - 1);
%endfunction
%}


 
%endfunction
  
%function Terminate(block)

%disp(['Terminating the block with handle ' num2str(block.BlockHandle) '.']);

%endfunction

是一个简单的模型。

enter description here

输出一个方波，时间间隔是0.01s，这里可以看出，block.Dwork(1).Data确实是一个全局变量。

回到Cassie电弧模型的离散化。是因为deltat这个变量的使用出现问题。

function msfunCassie_Rv_discrete(block)
% smileland作品，免费使用，但请不要删除此句话
% 以下设定参数并不一定全部需要，不需要的注释掉就好
setup(block);

%endfunction

function setup(block)
  %block.CurrentTime 只读对象，调用当前时间

  % 初始化对话框参数，即对话框内给定的初始化参数.
  block.NumDialogPrms     = 3;%给定3个参数
  %block.DialogPrmsTunable = {'Tunable','Nontunable','SimOnlyTunable'};%参数是否可变（此句话可选是否需要，默认不可变）

  %%%%%%%%初始化输入输出端口的属性%%%%%%%
  
  % 初始化输入输出端口个数.
  block.NumInputPorts  = 1;%输入端口个数
  block.NumOutputPorts = 1;%输出端口个数
  
  % 设定输入输出端口时继承还是动态，以指示输入和输出端口从模型继承其已编译属性（尺寸、数据类型、复杂度和采样模式）.
  block.SetPreCompInpPortInfoToDynamic;%设定输入端口为动态
  block.SetPreCompOutPortInfoToDynamic;%设定输出端口为动态

  % 设定输入端口属性.（也可以再Dwork中设定，在PostPropagationSetup子方法中）
  block.InputPort(1).DatatypeID  = 0;  % （-1 inherited；0 double;1 single;2 int8;3 uint8;4 int16;5 uint16;6 int32;7 uint12;8 boolean）
  block.InputPort(1).Complexity  = 'Real';%是否为复数
 
  % 设定输出端口属性.（也可以再Dwork中设定，在PostPropagationSetup子方法中）
  block.OutputPort(1).DatatypeID  = 0; % double
  block.OutputPort(1).Complexity  = 'Real';
  
  block.InputPort(1).Dimensions = 1;%输入口1的维数是1，最高支持到2维
  % block.InputPort(1).DimensionsMode = 'Fixed';%Fixed表示端口维数固定，Variable表示端口维数可变
  block.OutputPort(1).Dimensions = 1;%输出口1的维数是1，最高支持到2维
  % block.OutputPort(1).DimensionsMode = 'Fixed';%Fixed表示端口维数固定，Variable表示端口维数可变
  
  
  block.InputPort(1).DirectFeedthrough = true;%输入口有无直接馈通（true/false）

  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  % 设定有多少连续状态变量.
  %block.NumContStates = 1;
  
  % 设定有多少离散状态变量.
  % block.NumDworkDiscStates = 1;

  % 设定采样时间.
  %  [0 0.5]            : 从0.5s开始连续采样
  %  [0.1 0.5]          : 从0.5s开始，每0.1s采样一次
  %  [-1, 0]            : 继承采样频率
  %  [-2, 0]            : 动态变化的采样频率
  block.SampleTimes = [2.5e-5 0];%连续采样率
  
  %设定离散工作向量个数
  % block.NumDworks = 1;
  
  % -----------------------------------------------------------------
  % Options，若使用TLC则需要以下语句，否则注释掉即可
  % -----------------------------------------------------------------

  % block.SetAccelRunOnTLC(false);
  
  % block.OperatingPointCompliance = 'Default';
  
  % -----------------------------------------------------------------
  % 引用使用的子方法，主要有PostPropagationSetup;InitializeConditions;
  % Start;Outputs;Updates;Derivatives;Terminate.以下子方法仅保留着几个。
  % -----------------------------------------------------------------
  
  %
  %包含PostPropagationSetup子方法，用于初始化DWork工作向量
  block.RegBlockMethod('PostPropagationSetup', @DoPostPropSetup);
  %}
  
  %包含InitializeConditions子方法，用于初始化状态变量或DWork工作向量的值
  block.RegBlockMethod('InitializeConditions', @InitializeConditions);
  
  %与InitializeConditions子方法功能一致，但只初始化一次，Ini方法每次调用子方法时均会初始化
  % block.RegBlockMethod('Start', @Start);

  %包含输出模块
  block.RegBlockMethod('Outputs', @Outputs);

  %包含离散状态更新模块
  block.RegBlockMethod('Update', @Update);

  %包含连续状态微分模块
  %block.RegBlockMethod('Derivatives', @Derivatives);
  
  %包含终止模块
  %block.RegBlockMethod('Terminate', @Terminate);

  %还有其他模块，仅将最主要的列出，若需要则使用S函数自带模板。
%endfunction

% -------------------------------------------------------------------
% 以下为子方法模块，同样主要有DoPostPropSetup;InitializeConditions;
  % Start;Outputs;Updates;Derivatives;Terminate.以下子方法仅保留这几个。
% -------------------------------------------------------------------
%
function DoPostPropSetup(block)
  block.NumDworks = 4;
  
  block.Dwork(1).Name            = 'lngt1';
  block.Dwork(1).Dimensions      = 1;
  block.Dwork(1).DatatypeID      = 0;      % double
  block.Dwork(1).Complexity      = 'Real'; % real
  block.Dwork(1).UsedAsDiscState = true;
  
  block.Dwork(2).Name            = 't';
  block.Dwork(2).Dimensions      = 1;
  block.Dwork(2).DatatypeID      = 0;      % double
  block.Dwork(2).Complexity      = 'Real'; % real
  block.Dwork(2).UsedAsDiscState = true;
  
  block.Dwork(3).Name            = 'step';
  block.Dwork(3).Dimensions      = 1;
  block.Dwork(3).DatatypeID      = 0;      % double
  block.Dwork(3).Complexity      = 'Real'; % real
  block.Dwork(3).UsedAsDiscState = true;
  
  block.Dwork(4).Name            = 'lngt0';
  block.Dwork(4).Dimensions      = 1;
  block.Dwork(4).DatatypeID      = 0;      % double
  block.Dwork(4).Complexity      = 'Real'; % real
  block.Dwork(4).UsedAsDiscState = true;
  
  % Register all tunable parameters as runtime parameters.
  block.AutoRegRuntimePrms;

%endfunction
%}
%

function InitializeConditions(block)

block.Dwork(1).Data = log(1e-6);


    
%{
function Start(block)

  tau = block.DialogPrm(1).Data;
  Uc = block.DialogPrm(2).Data; 

%}

%
  function Update(block)

      tau = block.DialogPrm(1).Data;

      Uc = block.DialogPrm(2).Data; 
      
      if block.CurrentTime <= block.DialogPrm(3).Data
          step = 0;
      else
          step = 1;
          deltat = block.CurrentTime - block.Dwork(2).Data;
          block.Dwork(2).Data = block.CurrentTime;
      end
      
      block.Dwork(1).Data = step * 1/tau * (block.InputPort(1).data^2 / Uc^2 - 1) * 2.5e-5 + block.Dwork(4).Data;
      block.Dwork(4).Data = block.Dwork(1).Data;
      
%}

function Outputs(block)
  
  block.OutputPort(1).Data = block.InputPort(1).data * (exp(block.Dwork(1).Data));
  
%{
function Derivatives(block)

tau = block.DialogPrm(1).Data;
Uc = block.DialogPrm(2).Data; 

if block.CurrentTime <= block.DialogPrm(3).Data
    block.Dwork(3).Data = 0;
else
    block.Dwork(3).Data = 1;
end
    
block.Derivatives.Data = block.Dwork(3).Data * 1/tau * (block.InputPort(1).data^2 / Uc^2 - 1);
%endfunction
%}


 
%endfunction
  
%function Terminate(block)

%disp(['Terminating the block with handle ' num2str(block.BlockHandle) '.']);

%endfunction

一个不完全版的程序，可以实现2.5e-5的采样率的离散电弧模型。

enter description here

模型名称power_arcmodels_cassie_sfunction_discrete

现在可以回到单相接地故障模型中了。