求解输出调节方程

考虑如下系统
$$\{\begin{array}{l}\dot{x}=Ax+Bu\\ y=Cx\end{array}$$

$$\{\begin{array}{l}\dot{\upsilon }=S\upsilon \\ \xi =F\upsilon \end{array}$$

其中$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}$。

当我们对上述系统进行输出同步控制时(如异构多智能体的输出同步)，经常会使用到输出同步方程，即
$$\{\begin{array}{l}AX+BU=XS\\ CX=F\end{array}$$

上式中$X$和$U$便是我们所要求解的.
对于该方程直接求解可能会比较麻烦，我们可将其转化为线性方程的形式来进行求解.

考虑到更一般的形式，我们将上式转化为
$$\{\begin{array}{l}AX+BU-XS=E\\ CX=F\end{array}$$

这是因为在有些系统中会有$E$矩阵，若没有$E$的话，则可将$E$看作是零矩阵，并且其中各矩阵的维度为$E\in {\mathbb{R}}^{n\times q},X\in {\mathbb{R}}^{n\times q},U\in {\mathbb{R}}^{m\times q}$.

基于参考文献中的描述

我们可将输出同步方程写成如下的形式

$$\left(I_q\otimes \left[\begin{array}{cc}A& B\\ C& 0\end{array}\right]\right)\left[\begin{array}{c}{X}_1\\ U_1\\ ⋮\\ X_q\\ U_q\end{array}\right]-\left(S^T\otimes \left[\begin{array}{cc}{I}_n& 0\\ 0& 0\end{array}\right]\right)\left[\begin{array}{c}{X}_1\\ U_1\\ ⋮\\ X_q\\ U_q\end{array}\right]=\left[\begin{array}{c}{E}_1\\ F_1\\ ⋮\\ E_q\\ F_q\end{array}\right]$$

其中$\left[\begin{array}{cc}A& B\\ C& 0\end{array}\right]\in R^{(n+p)\times (n+m)}$,$X_q$表示$X$的第$q$列.

从而有$Q=I_q\otimes \left[\begin{array}{cc}A& B\\ C& 0\end{array}\right]-S^T\otimes \left[\begin{array}{cc}{I}_n& 0\\ 0& 0\end{array}\right]$（注意这边我推出来的和图片上的正负号是反的，这取决于具体的移项方法，是左移，还是右移，自己再计算的过程中，一定要注意）

进一步便可以得到输出同步方程转化为

$$Q\chi =b$$

其中$\chi ={\left[\begin{array}{ccccc}{X}_1^T& U_1^T& \cdots & X_q^T& U_q^T\end{array}\right]}^T$, $b={\left[\begin{array}{ccccc}{E}_1^T& F_1^T& \cdots & E_q^T& F_q^T\end{array}\right]}^T$.

取如下矩阵$A=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& -1& -2\end{array}\right]$, $B=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$, $C=\left[\begin{array}{ccc}0& 0& 1\end{array}\right]$, $S=\left[\begin{array}{cc}0& 1\\ -2& 0\end{array}\right]$, $F=\left[\begin{array}{cc}1& 0\end{array}\right]$, $E=\left[\begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}\right]$.

具体代码为

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];