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(:,j)
表示选出来的两列向量组成的矩阵,":"
表示行全选
我们要完成第二件事,其实只需要让 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
声明:未经本人同意,禁止转载、搬运、抄袭!
本文使用 markdown++进行格式转换
https://greasyfork.org/zh-CN/scripts/457903-bilibili-markdown