【数理知识】Riccati 黎卡提 system
相關鏈接:【Matlab】求解黎卡提 Riccati 方程 李雅普諾夫 Lyapunov 方程
黎卡提 Riccati
- 黎卡提方程 (Riccati equation)
- 代數 Riccati 方程
- 1 離散時間的代數Riccati方程
- 2 求解
- 離散代數黎卡提方程求解
- 1 黎卡提方程
- 2 迭代法求解
- 3 Matlab 仿真代碼
- 使用 ode45 求解 Riccati 方程
黎卡提方程 (Riccati equation)
黎卡提方程是最簡單的一類非線性方程。形如
y′=P(x)y2+Q(x)y+R(x)y'=P(x)y^2+Q(x)y+R(x)y′=P(x)y2+Q(x)y+R(x)
的方程稱為黎卡提方程。
From: 黎卡提方程-百度百科
代數 Riccati 方程
代數 Riccati 方程(algebraic Riccati equation)是最優控制的非線性方程,和連續時間或是離散時間下,無限時間(infinite-horizon)的最優控制有關。
標準的代數Riccati方程分為如下兩種:
連續時間代數Riccati方程(CARE):
ATP+PA?PBR?1BTP+Q=0A^TP + PA - PBR^{-1}B^TP + Q = 0ATP+PA?PBR?1BTP+Q=0
離散時間代數Riccati方程(DARE):
P=ATPA?(ATPB)(R+BTPB)?1(BTPA)+QP = A^TPA - (A^TPB)(R+B^TPB)^{-1}(B^TPA) + QP=ATPA?(ATPB)(R+BTPB)?1(BTPA)+Q
PPP 是未知數的 n×nn \times nn×n 對稱矩陣,A、B、QA、B、QA、B、Q 及 RRR 是已知實系數矩陣。
一般而言此方程式有許多的解,不過若有存在穩定解的話,希望可以找到穩定解。
1 離散時間的代數Riccati方程
在無限時間的最佳控制問題中,關注的是一些變數在相當時間之后的值,因此需在現在選定控制變數的值,讓系統在之后的時間都在最佳狀態下運作。控制變數在任意時間下的最佳值可以用 Riccati 方程的解以及狀態變數當時的觀測值求得。若觀測變數及控制變數都不只一個,Riccati 方程就會是矩陣方程。
代數 Riccati 方程可以決定無限時間下非時變 LQR 控制器的解,以及無限時間下非時變 LQG 控制的解。這兩個是控制理論中的基礎問題。
LQR 控制器
最優控制理論主要探討的是讓動力系統以在最小成本來運作,若系統動態可以用一組線性微分方程表示,而其成本為二次泛函,這類的問題稱為線性二次(LQ)問題。此類問題的解即為線性二次調節器,簡稱LQR。
LQG 控制
LQG控制 的全名是線性二次高斯控制,是控制理論中的基礎最優控制問題之一。此問題和存在加性高斯白噪聲的線性系統有關。此問題是要找到最佳的輸出回授律,可以讓二次費用函數的期望值最小化。其輸出量測假設受到高斯噪聲的影響,其初值也是高斯隨機向量。
控制理論
控制理論是工程學與數學的跨領域分支,主要處理在有輸入信號的動力系統的行為。系統的外部輸入稱為“參考值”,系統中的一個或多個變數需隨著參考值變化,控制器處理系統的輸入,使系統輸出得到預期的效果。
典型的離散時間 LQR 問題,是要最小化以下的函數:
∑t=1T(ytTQyt+utTRut)\sum_{t=1}^{T} (y_t^T Q y_t + u_t^T R u_t)t=1∑T?(ytT?Qyt?+utT?Rut?)
其狀態方程如下:
yt=Ayt?1+Buty_t = Ay_{t-1} + Bu_tyt?=Ayt?1?+But?
其中 yyy 是 n×1n × 1n×1 的狀態變數向量,uuu 是 k×1k × 1k×1 的控制變數向量,AAA 是 n×nn × nn×n 的狀態遞移矩陣,BBB 是 n×kn × kn×k 的控制系數矩陣,Q(n×n)Q (n × n)Q(n×n) 是對應半正定狀態損失函數矩陣,R(k×k)R (k × k)R(k×k) 是對應正定的控制損失函數矩陣。
從最后時間往前的推導可以找到每一個時間的最佳控制解:
ut?=?(BTPtB+R)?1(BTPtA)yt?1u_t^* = -(B^TP_tB + R)^{-1} (B^T P_t A) y_{t-1}ut??=?(BTPt?B+R)?1(BTPt?A)yt?1?
其中對應正定 cost-to-go(代價函數) 矩陣 PPP 會依據下式,配合 PT=QP_T = QPT?=Q,以逆向時間推導
Pt?1=Q+ATPtA?ATPtB(BTPtB+R)?1BTPtAP_{t-1} = Q + A^T P_t A - A^T P_t B (B^TP_t B + R)^{-1} B^T P_t APt?1?=Q+ATPt?A?ATPt?B(BTPt?B+R)?1BTPt?A
這個就是離散時間的代數 Riccati 方程。PPP 的穩態解和和 TTT 趨近無限大時的無限時間問題有關,可以將動態方程反復迭代直到收斂,來求得 PPP 的穩態解,之后再將動態方程中的時間標注移除,來確認穩態解是否正確。
2 求解
若代數 Riccati 方程存在穩定解,求解器一般會設法找到唯一的穩定解。穩定解的意思是指用此解控制相關的 LQR 系統,可以使閉回路的系統穩定。
針對CARE,其控制律為:
K=R?1BTPK = R^{-1} B^T PK=R?1BTP
閉回路遞移矩陣為:
A?BK=A?BR?1BTPA-BK = A-BR^{-1}B^TPA?BK=A?BR?1BTP
其穩定的充分必要條件是所有的特征值都有負的實部。
針對 DARE,其控制律為:
K=(R+BTPB)?1BTPAK = (R + B^TPB)^{-1} B^T PAK=(R+BTPB)?1BTPA
閉回路遞移矩陣為
A?BK=A?B(R+BTPB)?1BTPAA-BK = A - B(R+B^TPB)^{-1}B^TPAA?BK=A?B(R+BTPB)?1BTPA
其穩定的充分必要條件是所有的特征值在復數平面的單位圓內。
代數 Riccati 方程的解可以用 Riccati 方程的的迭代或是矩陣因式分解求得。離散時間問題的一種迭代方式是由有限時間問題下的動態 Riccati 方程,每一次迭代時,矩陣中的值都是從最終時間往前一段有限時間內的最佳解,若進行無限長的迭代。就會分斂到特定矩陣,是無限時間內的最佳解。
針對大型系統,也可以用找特征分解的方式求解。針對CARE,可以定義漢彌爾頓矩陣 (Hamiltonian):
Z=(A?BR?1BT?Q?AT)Z = \left(\begin{matrix} A & -BR^{-1}B^T \\ -Q & -A^T \end{matrix}\right)Z=(A?Q??BR?1BT?AT?)
因為 ZZZ 是漢彌爾頓矩陣 (Hamiltonian),若在虛軸上沒有特征值,則會有恰好一半的特征值會有負的實部。若定義 2n×n2n\times n2n×n 矩陣,其縱排 (column) 形成對應子空間的基底,表示為區塊矩陣的形式,如下所示:
(U1U2)\left(\begin{matrix} U_1 \\ U_2 \end{matrix}\right)(U1?U2??)
則
P=U2U1?1P = U_2 U_1^{-1}P=U2?U1?1?
是 Riccati 方程的解。而且 A?BR?1BTPA-BR^{-1}B^TPA?BR?1BTP 的特征值即為 ZZZ 特征值中有負實部的特征值。
針對 DARE,若 AAA 是可逆矩陣,可以定義辛矩陣:
Z=(A+BR?1BT(A?1)TQ?BR?1BT(A?1)T?(A?1)TQ(A?1)T)Z = \left(\begin{matrix} A+BR^{-1}B^T(A^{-1})^TQ & -BR^{-1}B^T(A^{-1})^T \\ -(A^{-1})^{T}Q & (A^{-1})^T \end{matrix}\right)Z=(A+BR?1BT(A?1)TQ?(A?1)TQ??BR?1BT(A?1)T(A?1)T?)
因為 ZZZ 是辛矩陣,若在單位圓圓周上沒有特征值,則會有恰好一半的特征值會在單位圓內。若定義 2n×n2n\times n2n×n 矩陣,其縱排 (column) 形成對應子空間的基底,表示為區塊矩陣的形式,如下所示:
(U1U2)\left(\begin{matrix} U_1 \\ U_2 \end{matrix}\right)(U1?U2??)
則
P=U2U1?1P = U_2 U_1^{-1}P=U2?U1?1?
是 Riccati 方程的解。而且 A?B(R+BTPB)?1BTPAA-B(R+B^TPB)^{-1}B^TPAA?B(R+BTPB)?1BTPA 的特征值即為 ZZZ 特征值中在單位圓內的特征值。
From: 代數Riccati方程
離散代數黎卡提方程求解
1 黎卡提方程
在LQR最優控制中,有連續時間最優控制,即LQR,也有離散時間最優控制DLQR,則在求解中一定會遇到解連續時間黎卡提方程和離散時間黎卡提方程的問題,本文主要針對離散時間黎卡提方程進行求解。
給定一個系統,如下所示:
x(k+1)=Ax(k)+Bu(k)x(k+1) = Ax(k)+Bu(k)x(k+1)=Ax(k)+Bu(k)
性能指標:
J(u)=∑n=1∞(xT(n)Qx(n)+uT(n)Ru(n)+2xT(n)Lu(n))J(u) = \sum_{n=1}^{\infty} (x^T(n) Q x(n) + u^T(n) R u(n) + 2x^T(n) L u(n))J(u)=n=1∑∞?(xT(n)Qx(n)+uT(n)Ru(n)+2xT(n)Lu(n))
離散時間黎卡提方程:
P=ATPA?ATPB(R+BTPB)?1BTPA+QK=(R+BTPB)?1BTPAu(n)=?Kx(n)P = A^TPA - A^TPB(R+B^TPB)^{-1}B^TPA + Q \\ K = (R+B^TPB)^{-1} B^TPA \\ u(n) = -K x(n)P=ATPA?ATPB(R+BTPB)?1BTPA+QK=(R+BTPB)?1BTPAu(n)=?Kx(n)
即求解上述黎卡提方程有兩種方法,一種是不變子空間法,另外一種是迭代法求解。
2 迭代法求解
P(k)=(A?BR?1LT)T(P?1(k?1)+BR?1BT)?1(A?BR?1LT)+Q?LR?1LTP(k) = (A-BR^{-1}L^T)^T (P^{-1}(k-1)+BR^{-1}B^T)^{-1} (A-BR^{-1}L^T) + Q - LR^{-1} L^TP(k)=(A?BR?1LT)T(P?1(k?1)+BR?1BT)?1(A?BR?1LT)+Q?LR?1LT
3 Matlab 仿真代碼
%======求解黎卡提方程==========% clear all;clc;close all; A = [0.8 0.3;-0.6 0]; B = [1 0.3;0 1.4]; L = [2 0.5;1 0.3]; Q = [7 1;1 3]; R = [3 1;1 2]; [Kr,Pr]=dlqr(A,B,Q,R,L); % 標準庫求解 err = 10^(-15); error = 10; Pe = Q - L*inv(R)*L'; X = Q - L*inv(R)*L'; % X = (A - B*inv(R)*L')'*inv(B*inv(R)*B')*(A-B*inv(R)*L') +Q - L*inv(R)*L'; Last_Pe = Pe; i = 1; while (error>=err)Pe = (A-B*inv(R)*L')'*inv(Last_Pe + B*inv(R)*B')*(A-B*inv(R)*L') +Q - L*inv(R)*L';error = norm((Pe-Last_Pe),'Inf');Last_Pe = Pe;i = i +1; end Ke = inv(R+B'*Pe*B)*B'*Pe*A;From: 求解離散黎卡提矩陣代數方程
使用 ode45 求解 Riccati 方程
ODE45——求解狀態變量(微分方程組)
總結
以上是生活随笔為你收集整理的【数理知识】Riccati 黎卡提 system的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 【控制】《多智能体系统的协同群集运动控制
- 下一篇: 【数理知识】拉格朗日乘数 Lagrang