MATLAB ode45耦合矩阵ODE

Question

我已经学习了最优控制中线性二次跟踪问题的两个耦合矩阵ODE，它们是：

​

哪里

​

​

我正在试着写一个能同时解微分方程的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正确解决耦合矩阵常微分方程的一般方法？

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的行重新构造矩阵。你可以用它来创建显式的助手函数。