解微分方程单时间循环利用MATLAB

Question

我有以下公式的机械系统：

xdot = Ax+ Bu

我要解决的一个循环这个等式，因为在我的每一步需要更新û但求解器如ode45或lsim解决时间间隔的微分方程。

for i = 1:10001 
    if x(i,:)>= Sin1 & x(i,:)<=Sout2 
     U(i,:) = Ueq - (K*(S/Alpha)) 
    else 
     U(i,:) = Ueq - (K*S) 
    end 
    % [y(i,:),t,x(i+1,:)]=lsim(sys,U(i,:),(time=i/1000),x(i,:)); 
    or %[t,x] = ode45(@(t,x)furuta(t,x,A,B,U),(time=i/1000),x) 
end 

难道我有另一种方法来解决这个方程在循环中的一次（不是单一的时间步）。

Answer 1

有许多方法可以跨函数调用更新和存储数据。 对于ODE套件，我已经开始喜欢所谓的“闭包”。 闭包基本上是一个嵌套函数，它从父函数中访问或修改一个变量。

下面的代码通过将传递给ode45和'OutputFcn'的右侧函数包装在名为odeClosure()的父函数中来使用此功能。

您会注意到我正在使用逻辑索引而不是if-声明。 if中的矢量 - 仅当所有元素都为真时，陈述才会成立，反之亦然。 因此，我创建一个逻辑数组并使用它来分母1或Alpha，具体取决于每行x/U的信号值。

'OutputFcn'storeU()在成功的时间步后ode45被调用。 该功能增长U存储阵列并适当更新它。 数组0123¾将具有与tspan（本例中为12）请求的解答点数相同的列数。 如果一个成功的完整步骤跳过了所有要求的点，该函数被调用时，中间所有请求的时间及其相关的解决方案值（因此x可能是矩形的，而不仅仅是一个矢量）。这就是为什么我使用bsxfun而不是rhs。

实例功能：

function [sol,U] = odeClosure() 

    % Initilize 
%  N  = 10   ; 
    A  = [ 0,0,1.0000,0; 0,0,0,1.0000;0,1.3975,-3.7330,-0.0010;0,21.0605,-6.4748,-0.0149]; 
    B  = [0;0;0.6199;1.0752 ] ; 
    x0 = [11;11;0;0]; 
    K  = 100; 
    S  = [-0.2930;4.5262;-0.5085;1.2232]; 
    Alpha = 0.2   ; 
    Ueq = [0;-25.0509;6.3149;-4.5085]; 
    U  = Ueq; 
    Sin1 = [-0.0172;-4.0974;-0.0517;-0.2993]; 
    Sout2 = [0.0172 ; 4.0974; 0.0517; 0.2993]; 

    % Solve 
    options = odeset('OutputFcn', @(t,x,flag) storeU(t,x,flag)); 
    sol  = ode45(@(t,x) rhs(t,x),[0,0.01:0.01:0.10,5],x0,options); 


    function xdot = rhs(~,x) 

     between = (x >= Sin1) & (x <= Sout2); 
     uwork = Ueq - K*S./(1 + (Alpha-1).*between); 
     xdot = A*x + B.*uwork; 

    end 

    function status = storeU(t,x,flag) 

     if isempty(flag) 
      % grow array 
      nAdd  = length(t)   ; 
      iCol  = size(U,2) + (1:nAdd); 
      U(:,iCol) = 0     ; 

      % update U 
      between = bsxfun(@ge,x,Sin1) & bsxfun(@le,x,Sout2); 
      U(:,iCol) = Ueq(:,ones(1,nAdd)) - K*S./(1 + (Alpha-1).*between); 
     end 

     status = 0; 
    end 

end