MATLAB ode45耦合矩阵ODE

Question

我已经了解了两个耦合矩阵ODE，它们用于最优控制中的线性二次跟踪问题，如下所示：

where

我正在尝试编写同时解决微分方程的MATLAB。这是我到目前为止的内容：

function [dSdt dGdt] = mRiccati2(t, S, A, B, Q, R, G, r, h)
    k = 1+floor(t/h);
    S = reshape(S, size(A)); %Convert from "n^2"-by-1 to "n"-by-"n"
    dSdt = A'*S + S*A - S*B*inv(R)*B'*S + Q; %Determine derivative
    dGdt = -(A'- S*B*inv(R)*B')*G + Q*r(:,k);
    dSdt = dSdt(:); %Convert from "n"-by-"n" to "n^2"-by-1
end

并且我尝试将函数调用为

[T S G] = ode45(@(t, S, G)mRiccati2(t, S, A, B, Q, R, G, r, h), [0:h:Tfinal], S0, G0);

不幸的是，我收到此错误：

Not enough input arguments.

Error in HW5 (line 24)
[T S G] = ode45(@(t, S, G)mRiccati2(t, S, A, B, Q, R, G, r, h), [0:h:Tfinal], S0, G0);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in HW5 (line 24)
[T S G] = ode45(@(t, S, G)mRiccati2(t, S, A, B, Q, R, G, r, h), [0:h:Tfinal], S0, G0);

是否有一般方法可以用ode45正确求解耦合矩阵ODE？

Answer 1

Matlab做了很多小事情，以使其更容易获得良好的结果，但这并不明智。您需要使问题适应Matlab的界面，没有自动检测功能。因此，您赋予ode45的函数需要消耗一个平面数组作为状态，并返回一个平面数组作为导数。您已经找到了进行转换的答案，只需要进行到底即可。

function dXdt = mRiccati2(t, X, A, B, Q, R, r, h)
    k = 1+floor(t/h);
    n = size(A(1,:))
    X = reshape(X, [n+1,n]); %Convert from flat to matrix, first the rows of S, then G
    S = X(1:n,:);
    G = X(n+1,:);
    dSdt = A'*S + S*A - S*B*inv(R)*B'*S + Q; %Determine derivative
    dGdt = -(A'- S*B*inv(R)*B')*G + Q*r(:,k);
    dXdt = [ dSdt(:) dGdt(:) ]; %Convert from matrix objects to flat arrays
end

然后当然要相应地调用积分器，根据初始矩阵将初始数据构造为平面数组

[T X] = ode45(@(t, X)mRiccati2(t, X, A, B, Q, R, r, h), [0:h:Tfinal], [S0(:) G0(:) ]);

要使用结果，您需要以与导数函数相同的方式从X的行中重构矩阵。您可以使用它来创建显式的辅助函数。