對于代數(shù)Riccati方程的求解網(wǎng)上能找到很多的資源，matlab也有成熟的函數(shù)，但是對于時(shí)變系統(tǒng)的Riccati矩陣微分方程，能找到的資料還比較少。

一、求解代數(shù)Riccati方程

可以在網(wǎng)上找到很多資料，如

• https://blog.csdn.net/m0_62299908/article/details/127807014

matlab也有相應(yīng)的一系列函數(shù) lqr、icare等。
對于這些函數(shù)不同的適用范圍自己目前了解的還不夠，之后補(bǔ)上。
這些函數(shù)到底能不能用于求解時(shí)變系統(tǒng)自己還沒搞清楚。

二、如何處理時(shí)變系統(tǒng)

參見matlab官方論壇 Solving Riccati differential equation with time variable matrix
這個(gè)帖子給出的方案是ode45中有一個(gè)關(guān)于處理時(shí)變項(xiàng)的例子。

matlab的 ode45 函數(shù)官方文檔的例子中有一個(gè) ODE with Time-Dependent Terms。
內(nèi)容如下

Consider the following ODE with time-dependent parameters
$$y^{\prime }(t)+f(t)y(t)=g(t).$$
The initial condition is $y_0=1$. The function f(t) is defined by the n-by-1 vector f evaluated at times ft. The function g(t) is defined by the m-by-1 vector g evaluated at times gt.

Create the vectors f and g.

ft = linspace(0,5,25);
f = ft.^2 - ft - 3;

gt = linspace(1,6,25);
g = 3*sin(gt-0.25);

Write a function named myode that interpolates f and g to obtain the value of the time-dependent terms at the specified time. Save the function in your current folder to run the rest of the example.
The myode function accepts extra input arguments to evaluate the ODE at each time step, but ode45 only uses the first two input arguments t and y.

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evaluate ODE at time t

Solve the equation over the time interval [1 5] using ode45. Specify the function using a function handle so that ode45 uses only the first two input arguments of myode. Also, loosen the error thresholds using odeset.

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

Plot the solution, |y|, as a function of the time points, |t|.

plot(t,y)

三、如何處理“矩陣”微分方程

在上一節(jié)中給出的方法實(shí)際上是已經(jīng)把矩陣微分方程變換成了一列，如果矩陣較大的話還是比較麻煩。
這個(gè)帖子給出了較好的解決方案 How can I solve the matrix Riccati differential equation within MATLAB?
或者還可參考
How can I solve a matrix differential equation within MATLAB?

給出的解決方案如下
The matrix Riccati differential equation:

 dX/dt = A'X + XA - XBB'X + Q

can be solved using the functions in the ODE suite.
Assume your dependent matrix “X” is “n”-by-“n”, and you have known “n”-by-“n” matrices “A”, “B”, and “Q”. The following method will solve the matrix Riccati differential equation. Save the following as a MATLAB file somewhere on the MATLAB Path.

function dXdt = mRiccati(t, X, A, B, Q)
X = reshape(X, size(A)); %Convert from "n^2"-by-1 to "n"-by-"n"
dXdt = A.'*X + X*A - X*B*B.'*X + Q; %Determine derivative
dXdt = dXdt(:); %Convert from "n"-by-"n" to "n^2"-by-1

Then, you can use the ODE45 function to solve this problem:

[T X] = ode45(@(t,X)mRiccati(t, X, A, B, Q), [0 10], X0) 

For example, using the sample data:

A = [1 1; 2 1]; 
B = [1; 1]; 
Q = [2 1; 1 1];
X0 = [1; 1; 1; 1];

You can use the following command to solve the system of differential equation

[T X] = ode45(@(t,X)mRiccati(t, X, A, B, Q), [0 10], X0) 

ODE45 returns “X” as a vector at each time step. You may use the following code to reshape each row of “X” to get the matrix and store it in a cell array:

[m n] = size(X);
XX = mat2cell(X, ones(m,1), n);
fh_reshape = @(x)reshape(x,size(A));
XX = cellfun(fh_reshape,XX,'UniformOutput',false);

The results of this can be verified by the LQR function:

[K,S,E] = lqr(A, B, Q, 1)

where “S” should have results very similar to the last elements in “X” or “XX”. The LQR function computes the steady-state value of the system. In this example, we generated the solution for up to “t = 10”, which is an adequate approximation of infinity for this problem.

For more information on ODE45 and other such solvers, refer to the function reference page for ODE45 in the MATLAB documentation.

對于該方法網(wǎng)上還有人追問，如有需要可查閱

• https://www.mathworks.com/matlabcentral/answers/225984-how-can-i-solve-matrix-riccati-differential-equation
• how can i solve matrix riccati differential equation?

四、關(guān)于逆向積分

在上面的兩個(gè)例子中給出的都是初始條件，但實(shí)際上Riccati方程給出的往往是終端條件，需要逆向(backward)積分。
對于簡單的情況，微分方程中并不含自變量。
此時(shí)逆向積分并不復(fù)雜，只需要在定義時(shí)間區(qū)間tspan時(shí)把初始時(shí)刻和終端時(shí)刻調(diào)換即可，參考鏈接為
https://www.mathworks.com/matlabcentral/answers/151336-solving-differential-riccati-equation-with-a-boundary-condition

示例如下

par = [1;2;1];  % q0, q1, and q2
yf = 1;
ti = 0; tf = 2;
opt = odeset('AbsTol',1.0e-07,'RelTol',1.0e-07);
[t,y] = ode45( @riccatiEquation, [tf,ti], yf ,opt, par);

% Visualize
plot(t,y)

function dydx = riccatiEquation(x,y,parameters)
q0 = parameters(1);
q1 = parameters(2);
q2 = parameters(3);

dydx = q0 + q1*y + q2*y*y;
end

或者采用平移的方式，可參考吳受章的最優(yōu)控制書。
如果實(shí)在不行，用變量代換轉(zhuǎn)換成初值問題是不是也行？沒有進(jìn)行仿真驗(yàn)證。

其實(shí)對于這類問題有一個(gè)更廣義的形式，即微分方程的多點(diǎn)邊值問題，對此該類問題不能使用ode45解決，而應(yīng)采用 BVP4C 或 BVP5C 來解決。參考鏈接如下

• Using ode45 to solve system of ODEs with some initial conditions and some terminal conditions
• Solve BVP with Multiple Boundary Conditions
• https://stackoverflow.com/questions/29162100/how-can-i-solve-a-system-of-odes-with-some-terminal-conditions-in-matlab

五、線性時(shí)變系統(tǒng)的Riccati矩陣微分方程

前段時(shí)間因事中斷了一段時(shí)間，現(xiàn)在自己已經(jīng)不面臨這個(gè)問題了…
按照前幾節(jié)的步驟似乎能解決這一問題，但沒有經(jīng)過仿真驗(yàn)證自己心里也不是很有底。

六、補(bǔ)充各類參考鏈接

時(shí)間有限，本來還想給幾個(gè)例子，看之后自己還有沒有機(jī)會吧。這里把一些自己在這個(gè)過程中的參考鏈接放上來

How can I solve riccati equation in Simulink with varying state space?

有一些鏈接給了代碼，沒來得及細(xì)看，但感覺似乎都是時(shí)不變系統(tǒng)的文章來源地址http://www.zghlxwxcb.cn/news/detail-751504.html

• 這個(gè)代碼聲稱解決的是時(shí)變矩陣微分方程，但得到的是近似解。而且還有論文支撐，看上去最靠譜 Recursive Approximate Solution to Matrix Diff. Riccati Eq.
• Solve Riccati Differential Equation (solve_riccati_ode)
• symmetric differential matrix Riccati equations
• 這份代碼和上一份是同一個(gè)人上傳到MathWorks網(wǎng)站上的，相差僅一天，可能是一樣的 The solve the matrix Riccati differential equation
• Algebraic Riccati Equation Solver
• 這份代碼也有論文，但似乎很專業(yè)，可能是數(shù)學(xué)專業(yè)的 Invariant Galerkin Ansatz Spaces and Davison-Maki Methods for the Numerical Solution of Differential Riccati Equations

到了這里，關(guān)于matlab求解時(shí)變系統(tǒng)的Riccati矩陣微分方程的文章就介紹完了。如果您還想了解更多內(nèi)容，請?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！