产品
LTI 即线性定常系统
以下只对能控性矩阵讲解
如果对于能控性矩阵不满秩，则可以进行能控性矩阵分解，把他变成能控部分和不能控部分
（能控也就指的是传递函数零极点不会对消）
matlab 提供了很多关于矩阵的函数，极大的方便了工程运算
比如 tf、ss、zpk 空间模型快速转换、rank 求秩、rref 求最大线性无关组，今天就用 matlab 提供的函数来实现一下能控能观系统分解
首先写出一个求能控性矩阵的 Uc 的判断
输入参数是 A,B,C,D
n=length(A);
if length(A)==1
disp("你认真的？")
elseif length(A)==2
Uc=[B A*B];
Uo=[C;C*A];
elseif length(A)==3
Uc=[B A*B (A^2)*B];
Uo=[C;C*A;C*(A^2)];
elseif length(A)==4
Uc=[B A*B (A^2)*B (A^3)*B];
Uo=[C;C*A;C*(A^2) C*(A^3)];
end
求好了 Uc 矩阵了可以用
rank(Uc)求出秩 和 length(A)做比较
不满秩就可以进行下一步分解了
假设这里是三维矩阵
我们知道要分解首先要找到这么一个向量，即从 Uc 中选取两列（维数-1）线性无关向量，再自己随意补上一个与他们俩都线性无关的列向量组成 Tc 矩阵
这里的难点就是如何找出两列线性无关的向量，并确保补上去的向量满足要求 我们当然可以自己写一个化简阶梯来判断，但 matlab 提供了一个 rref 函数
I=eye(n);
[~,j]=rref(Uc); %先选出Uc本身线性无关的维数的n-1列
[~,j2]=rref([Uc(:,j),I]); %再选出与这线性无关的最后一列返回的 j 参数就表示 Uc 极大线性无关组所在的的列数的最小值 这样我们就可以轻松完成第一件事，从 Uc 选出两列线性无关的向量啦
这里的表示选出来的两列向量组成的矩阵，表示行全选
我们要完成第二件事，其实只需要让 Uc 和单位阵 I 组成增广矩阵来找其中的极大线性无关组 还记不记得之前说过 rref 只会返回尽可能列数小的 j，这就可以保证返回的 j2 表示的列数，一定会有原来两列 Uc 和某一列单位阵，这样我们对增广矩阵进行j2列索引就可以拿到变换矩阵Tc啦
Augmentation=[Uc(:,j),I];
Tc= Augmentation(:,j2);拿到Tc后就很简单了，只需要对n-1的行和列进行索引就可以找到可控的子系统啦
其中A2，B2，C2代表变换后的空间表达式参数
A3，B3，C3代表能控子系统的空间表达式参数
Tc_inv=inv(Tc);
A2=Tc_inv*A*Tc;
B2=Tc_inv*B;
C2=C*Tc;
A3=A2(1:RankUc,1:RankUc);
B3=B2(1:RankUc,:);
C3=C2(:,1:RankUc);对于能观性矩阵的分解，只需要把矩阵的转置传入rref中就可以啦！由于大体相同，不再赘述
以下是完整代码，请到2019以上的matlab版本，实时脚本中运行
clear;clc
disp("请输入你的状态空间参数A,B,C,D，例如")
disp("A=[-2 2 -1;0 -2 0;1 -4 0];B=[0;0;1];C=[1,-1,1];D=[0];")
A=[-2 2 -1;
0 -2 0;
1 -4 0];
B=[ 0;
0;
1];
C=[1 -1 1];
D=0
[Tc,Tc_inv,Ac2,Bc2,Cc2,Ac3,Bc3,Cc3,Uc,RankUc]=deComConMa(A,B,C,D); %控制分解
Tc_inv
[To,To_inv,Ao2,Bo2,Co2,Ao3,Bo3,Co3,Uo,RankUo]=deComObsMa(A,B,C,D); %观测分解
To_inv
disp("注：变换后的状态子空间由于选取不同，有无数种可能，但是子系统一定是唯一的")
function [Tc,Tc_inv,A2,B2,C2,A3,B3,C3,Uc,RankUc]=deComConMa(A,B,C,D)
n=length(A);
if length(A)==1
disp("你认真的？")
elseif length(A)==2
Uc=[B A*B];
% Uo=[C;C*A];
elseif length(A)==3
Uc=[B A*B (A^2)*B];
% Uo=[C;C*A;C*(A^2)];
elseif length(A)==4
Uc=[B A*B (A^2)*B (A^3)*B];
% Uo=[C;C*A;C*(A^2) C*(A^3)];
end
RankUc=rank(Uc);
if RankUc<n
disp("Uc不满秩,能进行分解");
I=eye(n);
[~,j]=rref(Uc); %先选出Uc本身线性无关的维数的n-1列
[~,j2]=rref([Uc(:,j),I]); %再选出与这线性无关的最后一列
Augmentation=[Uc(:,j),I];
Augmentation(:,j2);
Tc= Augmentation(:,j2);
Tc_inv=inv(Tc);
A2=Tc_inv*A*Tc;
B2=Tc_inv*B;
C2=C*Tc;
A3=A2(1:RankUc,1:RankUc);
B3=B2(1:RankUc,:);
C3=C2(:,1:RankUc);
disp("Uc的秩为:");
disp(RankUc);
disp("Uc:");
disp(Uc);
disp("Tc:");
disp(Tc);
disp("变换后的状态空间表达式A2,B2,C2分别为:");
disp(A2);
disp(B2);
disp(C2);
disp("可控子系统系数A3,B3,C3分别为:")
disp(A3)
disp(B3)
disp(C3)
else
disp("Uc满秩,能控,无需分解")
end
end
function [To,To_inv,A2,B2,C2,A3,B3,C3,Uo,RankUo]=deComObsMa(A,B,C,D)
n=length(A);
if length(A)==1
disp("你认真的？")
elseif length(A)==2
% Uc=[B A*B];
Uo=[C;C*A];
elseif length(A)==3
% Uc=[B A*B (A^2)*B];
Uo=[C;C*A;C*(A^2)];
elseif length(A)==4
% Uc=[B A*B (A^2)*B (A^3)*B];
Uo=[C;C*A;C*(A^2) C*(A^3)];
end
RankUo=rank(Uo);
if rank(Uo)<n
disp("Uo不满秩,能进行分解");
I=eye(n);
[~,j]=rref(Uo');
[~,j2]=rref([Uo(j,:);I]');%再选出与这线性无关的最后一列
Augmentation=[Uo(j,:);I];
To_inv=Augmentation(j2,:);
To=inv(To_inv);
A2=To_inv*A*To;
B2=To_inv*B;
C2=C*To;
A3=A2(1:RankUo,1:RankUo);
B3=B2(1:RankUo,:);
C3=C2(:,1:RankUo);
disp("Uo的秩为:")
disp(RankUo)
disp("Uo:")
disp(Uo)
disp("To:")
disp(To)
% disp("To的逆:")
% disp(To_inv)
disp("变换后的状态空间表达式A2,B2,C2分别为:")
disp(A2)
disp(B2)
disp(C2)
disp("能观子系统系数A3,B3,C3分别为:")
disp(A3)
disp(B3)
disp(C3)
else
disp("Uo满秩,能观,无需分解");
end
end
作者：Onion
邮箱：bigonion@bigonion.cn
声明：未经本人同意，禁止转载、搬运、抄袭！
NameSpace：https://bigonion.cn
Origin:https://bigonion.cn/blog
本文使用markdown++进行格式转换
https://greasyfork.org/zh-CN/scripts/457903-bilibili-markdown
155-2731-8020
