首页<<

N维LTI系统的能控、能观性分解的matlab实现


前言

matlab 提供了很多关于矩阵的函数,极大的方便了工程运算
比如 tf、ss、zpk 空间模型快速转换、rank 求秩、rref 求最大线性无关组,今天就用 matlab 提供的函数来实现一下能控能观系统分解


LTI 即线性定常系统

以下只对能控性矩阵讲解

如果对于能控性矩阵不满秩,则可以进行能控性矩阵分解,把他变成能控部分和不能控部分

(能控也就指的是传递函数零极点不会对消)

首先写出一个求能控性矩阵的 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

2023-01-28 01:50

作者:Onion

邮箱:bigonion@bigonion.cn

声明:未经本人同意,禁止转载、搬运、抄袭!

本文使用 markdown++进行格式转换
https://greasyfork.org/zh-CN/scripts/457903-bilibili-markdown