考虑如下系统
{
x
˙
=
A
x
+
B
u
y
=
C
x
\begin{cases} \dot{x} =Ax+Bu\\ y=Cx \end{cases}
{x˙=Ax+Buy=Cx
{ υ ˙ = S υ ξ = F υ \begin{cases} \dot{\upsilon } =S\upsilon \\ \xi =F\upsilon \end{cases} {υ˙=Sυξ=Fυ
其中 A ∈ R n × n , B ∈ R n × m , C ∈ R p × n , S ∈ R q × q , F ∈ R p × q A\in \mathbb{R}^{n\times n} ,B\in \mathbb{R}^{n\times m} ,C\in \mathbb{R}^{p\times n} ,S\in \mathbb{R}^{q\times q} ,F\in \mathbb{R}^{p\times q} A∈Rn×n,B∈Rn×m,C∈Rp×n,S∈Rq×q,F∈Rp×q。
当我们对上述系统进行输出同步控制时(如异构多智能体的输出同步),经常会使用到输出同步方程,即
{
A
X
+
B
U
=
X
S
C
X
=
F
\begin{cases} AX+BU=XS\\ CX=F \end{cases}
{AX+BU=XSCX=F
上式中
X
X
X和
U
U
U便是我们所要求解的.
对于该方程直接求解可能会比较麻烦,我们可将其转化为线性方程的形式来进行求解.
考虑到更一般的形式,我们将上式转化为
{
A
X
+
B
U
−
X
S
=
E
C
X
=
F
\begin{cases} AX+BU-XS=E\\ CX=F \end{cases}
{AX+BU−XS=ECX=F
这是因为在有些系统中会有 E E E矩阵,若没有 E E E的话,则可将 E E E看作是零矩阵,并且其中各矩阵的维度为 E ∈ R n × q , X ∈ R n × q , U ∈ R m × q E\in \mathbb{R}^{n\times q} ,X\in \mathbb{R}^{n\times q} ,U\in \mathbb{R}^{m\times q} E∈Rn×q,X∈Rn×q,U∈Rm×q.
基于参考文献中的描述
我们可将输出同步方程写成如下的形式
( I q ⊗ [ A B C 0 ] ) [ X 1 U 1 ⋮ X q U q ] − ( S T ⊗ [ I n 0 0 0 ] ) [ X 1 U 1 ⋮ X q U q ] = [ E 1 F 1 ⋮ E q F q ] \left( I_{q} \otimes \begin{bmatrix} A & B\\ C & 0 \end{bmatrix}\right)\begin{bmatrix} X_{1}\\ U_{1}\\ \vdots \\ X_{q}\\ U_{q} \end{bmatrix} -\left( S^{T} \otimes \begin{bmatrix} I_{n} & 0\\ 0 & 0 \end{bmatrix}\right)\begin{bmatrix} X_{1}\\ U_{1}\\ \vdots \\ X_{q}\\ U_{q} \end{bmatrix} =\begin{bmatrix} E_{1}\\ F_{1}\\ \vdots \\ E_{q}\\ F_{q} \end{bmatrix} (Iq⊗[ACB0])⎣⎢⎢⎢⎢⎢⎡X1U1⋮XqUq⎦⎥⎥⎥⎥⎥⎤−(ST⊗[In000])⎣⎢⎢⎢⎢⎢⎡X1U1⋮XqUq⎦⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎡E1F1⋮EqFq⎦⎥⎥⎥⎥⎥⎤
其中 [ A B C 0 ] ∈ R ( n + p ) × ( n + m ) \begin{bmatrix} A & B\\ C & 0 \end{bmatrix} \in R^{( n+p) \times ( n+m)} [ACB0]∈R(n+p)×(n+m), X q X_{q} Xq表示 X X X的第 q q q列.
从而有 Q = I q ⊗ [ A B C 0 ] − S T ⊗ [ I n 0 0 0 ] Q=I_{q} \otimes \begin{bmatrix} A & B\\ C & 0 \end{bmatrix} -S^{T} \otimes \begin{bmatrix} I_{n} & 0\\ 0 & 0 \end{bmatrix} Q=Iq⊗[ACB0]−ST⊗[In000](注意这边我推出来的和图片上的正负号是反的,这取决于具体的移项方法,是左移,还是右移,自己再计算的过程中,一定要注意)
进一步便可以得到输出同步方程转化为
Q χ = b Q\chi =b Qχ=b
其中 χ = [ X 1 T U 1 T ⋯ X q T U q T ] T \chi =\begin{bmatrix} X^{T}_{1} & U^{T}_{1} & \cdots & X^{T}_{q} & U^{T}_{q} \end{bmatrix}^{T} χ=[X1TU1T⋯XqTUqT]T, b = [ E 1 T F 1 T ⋯ E q T F q T ] T b=\begin{bmatrix} E^{T}_{1} & F^{T}_{1} & \cdots & E^{T}_{q} & F^{T}_{q} \end{bmatrix}^{T} b=[E1TF1T⋯EqTFqT]T.
取如下矩阵 A = [ 0 1 0 0 0 1 0 − 1 − 2 ] A=\begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ 0 & -1 & -2 \end{bmatrix} A=⎣⎡00010−101−2⎦⎤, B = [ 0 0 1 ] B=\begin{bmatrix} 0\\ 0\\ 1 \end{bmatrix} B=⎣⎡001⎦⎤, C = [ 0 0 1 ] C=\begin{bmatrix} 0 & 0 & 1 \end{bmatrix} C=[001], S = [ 0 1 − 2 0 ] S=\begin{bmatrix} 0 & 1\\ -2 & 0 \end{bmatrix} S=[0−210], F = [ 1 0 ] F=\begin{bmatrix} 1 & 0 \end{bmatrix} F=[10], E = [ 0 0 0 0 0 0 ] E=\begin{bmatrix} 0 & 0\\ 0 & 0\\ 0 & 0 \end{bmatrix} E=⎣⎡000000⎦⎤.
具体代码为
A = [0 1 0; 0 0 1; 0 -1 -2];
B = [0; 0; 1];
C = [0 0 1];
S = [0 1; -2 0];
F = [1 0];
E = [0 0; 0 0; 0 0];
% 获取各矩阵的维度信息
[n, n] = size(A);
[q, q] = size(S);
[n, m] = size(B);
[p, n] = size(C);
% 线性方程中的矩阵
Q = kron(eye(q),[A B; C zeros(p,m)]) - kron(S', [eye(n), zeros(n,m); zeros(p,n) zeros(p,m)]);
b = [E(:,1);F(:,1);E(:,2);F(:,2)]; % 这里时因为选择的S只有两维
xx = linsolve(Q, b);
X = [xx(1:3) xx(5:7)]; % 1:3共3维,是因为X矩阵是2X3的
U = [xx(4) xx(8)]; % 4共1维,是因为U矩阵是1X2的
% 检验
X*S-A*X-B*U
C*X-F
计算出的解为
X=[-0.5 0; 0 -0.5; 1 0]
U=[2 0.5]
这里需要注意X,U的取值需要根据自己当前系统所得的维护来具体确定
本人只是根据文献中的阅读提供一个思路,上述代码直接复制可能会不适用,若有感兴趣,可写一个更通用的程序 ^_~
下面提供一个博主自己在进行仿真时,考虑不同的维度,用简单的if…else写的代码
function X_U = output_regulation(A, B, E, C, F, S, dim)
[~,m]=size(B);
[p,n]=size(C);
[~,q]=size(S);
Q = kron(eye(q),[A B; C zeros(p,m)]) - kron(S',[eye(n) zeros(n,m);zeros(p,n) zeros(p,m)]);
if dim == 2
b = [E(:,1);F(:,1);E(:,2);F(:,2)];
xx = linsolve(Q, b);
X = [xx(1:2) xx(4:5)];
U = [xx(3) xx(6)];
elseif dim == 3
b = [E(:,1);F(:,1);E(:,2);F(:,2)];
xx = linsolve(Q, b);
X = [xx(1:3) xx(5:7)];
U = [xx(4) xx(8)];
end
X_U = [X; U];