解微分方程單時間循環利用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