状态估计用于描述符 LTI 和 LPV 系统的分析状态估计和故障检测的算法（Matlab代码实现）

Posted 2023-04-02 我爱Matlab编程

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了状态估计用于描述符 LTI 和 LPV 系统的分析状态估计和故障检测的算法（Matlab代码实现）相关的知识，希望对你有一定的参考价值。

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码、数据、文章讲解

💥1 概述

文献来源：

摘要:本文讨论了在LTI或LPV描述子框架中建模的系统的MATLAB/SCILAB工具箱的开发。给出了正则性、可解性、可控性和可观测性。包括全阶观测器和降阶观测器，比例观测器和比例积分观测器。其中一些观察人士考虑了未知的输入。可以作为基于状态估计和故障检测的观测器的辅助工具。这些观察者已经从最近发表的几篇论文中得到了考虑。通过建立观测库来实现故障的检测和隔离。这些观察者库可以通过选择输入/输出矩阵来建立，也可以使用所提出的算法自动建立。

许多系统都可以用非线性微分方程来建模，但是监控系统的设计是一项困难的任务。因此，对非线性系统进行线性化以获得线性时不变(LTI)系统是非常常见的，但这种表示在一个平衡点或工作点附近是有效的。

改进模型表示的一种方法是包含一些限制。如果限制是模型的一部分，那么系统就变成了描述符- lti (DLTI)表示。DLTI系统的主要优点是集成了静态关系(例如物理约束)和动态关系。这些考虑允许对广泛的过程进行建模，例如，在Dai, 1989;段,2010)。

Luenberger观测者存在的充分条件在(Hou and Muller, 1999;Darouach和Boutayeb, 1995)。(Darouach等人，1996)的作者提出了一种降阶未知输入观测，类似于(Chen和Patton, 1999)中研究的LTI系统的观测器。

详细文章讲解及数学模型见第4部分。

📚2 运行结果

部分代码：

xe=[ 1 5 3 0]'% intial stimated states
ye(:,1)=C*xe % compute the initial output
t(1)=0;% initial time
u(1)=0;

% observer gains
% % if the system is observable then compute the restricted system equivalence systems by QR descomposition
disp('Simulation')

for k=1:50/Te
% Time counter
t(k+1)=t(k)+Te;
% input
if t(k)<7
u(k+1)=1;
else
u(k+1)=5;
end
% % the differentials equations can be solved by runge-kuta of four order
x1=x1+Te*(x3);
x2=x2+Te*(x1);
x4=x1;
x3=-x2-x4+u(k);
xx(:,k+1)=[x1 x2 x3 x4]';
y(:,k+1)=C*xx(:,k+1); % system output
% Observer
yi=[-Be1*u(k);y(:,k)]; % auxiliar for compute the observer
% Reduced order observer
% $z=Piz+Lyi+H*u $
% $\\hatx=Mz+Fyi$
z(k+1)=z(k)+Te*(Pi*z(k)+L*yi+H*u(k)); %
xe(:,k+1)=M*z(k)+F*yi;

ye(:,k+1)=C*xe(:,k+1); % stimated output

% Generation of residuos
end

disp('PLOTS')
pause
figure(2);
cla();
plot(t,y',t,ye')
legend('output', 'stimated')

n=size(A,1);
m=size(B1,2);
p=size(C,1);
if flag==1
P=sdpvar(n,n);
Y=sdpvar(m,n,'full');
Q=sdpvar(1,n,'full');
sdpvar g
% g=0.1
alpha=0;
Con=[P>=0,g>=0];
LMI=blkvar();
LMI(1,1)= (P*E'+S*Q)'*A'+Y'*B1'+A*(P*E'+S*Q)+B1*Y+B*B'+2*alpha*P;
LMI(1,2)=E*P'*C'+Q'*S'*C'+Y'*B2';
% (P*E'+S*Q)'*C'+Y'*B2';
LMI(2,2)=-g*eye(p);
LMI=sdpvar(LMI);
Con=[Con, LMI<=0]
op=sdpsettings('verbose',0,'solver','sedumi','sedumi.eps',1e-5);
solvesdp(Con,[],op)

P=double(P);
Y=double(Y);
Q=double(Q);
g=double(g);
% K=Y/((P*E'+S*Q))
K=Y/((P*E'+S*Q))
lamda=deig(A+B1*K,E)
% [Ccon, Rcon, Icon]=dcontr (E,A+B*K,B)
elseif flag==2
% g=0.01100
cvx_begin sdp
cvx_precision default
cvx_solver sedumi
variable g
variable P(n,n) symmetric
variable Y(m,n)
variable Q(1,n)
minimize g
P >= 0
[ (P*E'+S*Q)'*A'+Y'*B1'+A*(P*E'+S*Q)+B1*Y+B*B' E*P'*C'+Q'*S'*C'+Y'*B2' ; ...
(E*P'*C'+Q'*S'*C'+Y'*B2')' -g*eye(p)] <= 0
cvx_end

gamma = sqrt(g)
K=Y/((P*E'+S*Q));
lamda=deig(A+B1*K,E)

else

setlmis([]);
% create a blank LTI framework
P=lmivar(1,[n 1]); % declare X as a 3 脳 3 symmetrical matrix
Y=lmivar(2,[m n]); % declare Y as a n x n
Q=lmivar(2,[1 n]); % declare
lmiterm([1 1 1 P],A,E','s') % A*P*E' + *
lmiterm([1 1 1 Q],A*S,1,'s') % A*S*Q+*
lmiterm([1 1 1 Y],B1,1,'s') % B*Y+*
lmiterm([1 1 1 0],B*B') % B*B'+*
lmiterm([1 1 2 P'],E,C') % E*P'*C'
lmiterm([1 1 2 Q'],1,S'*C') % Q'*S'*C'
lmiterm([1 1 2 Y'],1,B2') % E*P'*C'
% lmiterm([1 1 1 P],2*5,1,'s') % B*Y+*
lmiterm([-2 1 1 P],1,1,'s') % P>0
% lmiterm([-3 1 1 Q],1,1) % P>0
Con=getlmis;
% c = mat2dec(Con,eye(3),eye(3),eye(1,n));
% options = [1e5,0,0,0,0];
% [copt,b] = mincx(Con,c,options);
[tmin b]=feasp(Con);
% [tmin b]=mincx(Con,-trac);
P=dec2mat(Con,b,P);
Y=dec2mat(Con,b,Y);
Q=dec2mat(Con,b,Q);
K=Y/((P*E'+S*Q));
lamda=deig(A+B*K,E);

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

🌈4 Matlab代码、数据、文章讲解

数字信号入门笔记2 —线性时不变（LTI）系统

数字信号入门笔记2 —线性时不变（LTI）系统

2.1系统与线性时不变系统

在信号处理背景下，系统可以定位为对输入信号进行某种处理，实现某种功能的物理结构。若用x(n)表示系统的输入，y(n)表示系统的输出，T[.]表示系统。对x(n)施加处理，得到y(n)，输入和输出可以表示为：

y(n) = T[x(n)]

一个系统任意时刻的输出取决于本时刻的输入，不依赖过去和将来时刻的输入，成为静态系统。放大器是一个静态系统。其他情况下，系统称为动态系统或有记忆系统，比如单位延时器。

2.1.1线性系统

对于一个系统T[.]，若输入为x1(n)时对应输出为y1(n)，输入为x2(n)时对应输出为y2(n)，即

y1(n)=T[x1(n)]，y2(n)=T[x2(n)]

满足齐次性: T[ax1(n)] = aT[x1(n)]

满足可加性: T[x1(n) + x2(n)] = T[x1(n)] + T[x2(n)]

上述两式可以合并为: T[ax1(n) + bx2(n)] = aT[x1(n)] + bT[x2(n)]，如果系统满足这个式子，则为线性系统。

2.1.2时不变系统

系统输入与输出关系不随时间变化而变化，或者说系统对输入信号的响应与加于系统的时间无关，则为时不变系统。

y(n-k) = T[x(n-k)]

实际中，可以先输入x(n)，计算出输出y(n)，然后将x(n)延时k个单位时间，再计算系统输出，如果两次输出序列相等，则为时不变系统。

2.2LTI系统的时域描述

2.2.1差分方程表示

这里我是这样理解，当前时刻的输出，既与当前n时刻和n时刻之前的输入有关，也与之前的输出有关。这里注意联系后文还会提到信号与系统建的作用方式为卷积，卷积会累加之前时刻的输入与系统间的作用值。

2.2.2单位冲击响应表示

单位冲击响应是输入信号为单位冲击信号ζ(n)时对应的系统输出，常用h(n)来表示。为什么LTI系统可以用单位冲击响应表示？单位冲击响应是一类最简单的信号，我们希望把一个复杂的信号分为多个类似ζ(n)的简单信号，通过考察系统对ζ(n)这样简单信号的输出来认识系统的特性。下来分两步来推导一下单位冲击信号如何表示系统。

1.我们先将一个离散信号x(n)分解为多个单位冲击信号之和

x(k) = x(k)ζ(n-k) = x(n)ζ(n-k)

这样系统在n=k时取值为x(k)，其他均为0，这样反复操作，就得到了x(n)。

ζ(n)只有在ζ(0)的时候才有值，所以x(n)需要和ζ(n-k)即ζ(0)相乘，得到在k个点时信号用单位冲击信号表示为 x(n)ζ(n-k)。这样就把一个离散信号信号分解为多个单位冲击信号相加：

2.将单位冲击信号通过一个LTI系统T[.]

根据LTI系统的齐次性，相当于系统先与单位冲击信号作用，得到一个系统的单位冲击响应，既h(n-m)。这里相当于x(n)与系统的单位冲击响应h(n-m)做卷积，最终得到输出y(n)。这样就由单位多个单位冲击信号叠加得到一个信号，将这个信号通过一个LTI系统，由LTI系统的齐次性，过度到了求出一个系统的单位冲击响应，信号与这个系统的单位冲击响应卷积得到输出y(n)，得出结论h(n)就是系统的结论。同时h(n)本身代表系统的单位冲击响应，这样就讲系统与系统的单位冲击相应联系在了一起，他们是等效的。后文还会提到FIR和IIR两种系统（滤波器），他们分别是单位冲击响应长度有限和单位冲击响应长度无限的系统。

2.3 LTI系统的特征信号复正弦信号

接上一节，已经得到了LTI系统的单位冲击响应表达式，如果对其施加一个复正弦信号，结果如下：

若k=n-m

可见这里输入的复正弦信号可以单独提出来，而其他部分变成了一个与x(n)无关的常数项（参考指数积分公式 ∫e^x dx = e^x+c）。换句话说就是对于不同角频率的复正弦信号，LTI系统只相当于其乘以一个常数（对于不同的角频率w0的信号这个常数不同），而输出信号的频率不变。由上推广，对于复杂的信号，可以讲复杂信号分解为多个复正弦信号，对于不同频率的复正弦信号，系统会输出不同的频率响应，这些频率响应就可以来表征一个系统。反过来看一下单位冲击信号的频谱：

这样，对一个系统输入一个单位冲击响应，在频域看相当于对其输入了所有频率的信号而得到的响应。这样看时域上信号分解为多个单位冲击响应和频域上信号分解为多个复正弦信号，是等效的两种操作。

2.4 Z变换分析LTI系统

1.Z变换的作用是方便的看出一个系统的性质，通过z变换可以在零极点图上描述系统，方便地看出系统是否稳定（单位冲击响应是否收敛)

2.Z变换提供了一种新的表示系统的方式，传递函数，只需要知道输入与输出序列的Z变换，就可以得到系统的传递函数

3.当z=e^jw 也就是在单位圆上时，这时的z变换就相当于傅里叶变换，这就是Z域与频域的关系。

2.4.1 Z变换定义

这里X(z)是一个无穷级数求和，x(n)的Z变换不一定总是存在。X(z)有限时，z的取值范围成为收敛域。下面是一些典型信号的Z变换，利用下表可以典型信号的Z变换可以方便推出常见x(n)的Z变换。

Z变换的反变换

这里记得大学复变函数考试基本上都是做这种复平面上封闭曲线积分运算，但是完全不知道这个东西有什么用，这里可以看出这种计算可以把一个Z变换后的结果反推出原始的时域信号。

2.4.2传递函数

这里需要先记住一个结论，两个序列时域的卷积在Z域会变为乘积。下面是对一个系统的输出进行Z变换，以此来引出传递函数。

展开Z变换

两个自变量变化相互独立，可以将积分拆开，变为两个积分相乘

将k=n-m

最终输出的Z变换Y(z)等于输入x(n)的Z变换X(z)与系统h(n)的Z变换H(z)的积。这样只要知道了系统的输入和输出，就可以得到系统的Z变换H(z)，经过Z反变换就可以得到h(n)，也就是系统的单位冲击响应。

2.4.3 Z变换单位延迟与差分方程

上一节通过传递函数和Z反变换建立了Z变换与系统单位冲击响应的关系，那么Z变换和差分方程间的关系是什么呢？

Z变换有一个移位特性

这个式子就建立了Z变换中第n个输出与第n-1次输出的关系。z^-1可以表示一个单位延时，有了单位延时就可以建立两次输出之间的关系，而差分方程同样也是代表前后两次间的关系。扩展到k个单位延时

这样现在将差分方程y(n)两边做Z变换

差分方程的Z变换

这样通过系统的差分方程就可以得到系统传递函数H(z)，进而得到h(n)。差分方程——(Z变换的单位延时)——>系统传递函数H(z)—(Z反变换)-—>系统单位冲击响应h(n) 三者都可以表征一个系统

2.4.4零极图直观体现系统特性

得到系统传递函数H(z)后，可以将其改写为如下形式是

满足N(z)等于0的zk，称为系统的零点，满足D(z)等于0的zk称为系统的极点。系统除了常数K，其他都是由零点极点决定的。由于系统函数H(z)的分子分母各项系数均为实数，所以零极点若为虚数或复数，则一定共轭成对出现。(这个地方没有看懂)

从零极图看单位冲击响应

有这样一个有一个实极点的系统

由Z变换表可知这个系统的单位冲击响应为

这是一个指数函数和一个单位阶跃函数的积。随着a的值的变化，系统的单位冲击响应的收敛程度是不一样的，具体如下。极点实数的情况下

下面看一下极点为复数的情况

对应单位冲击响应为

h(n)的形状主要由极点决定，零点主要影响h(n)的幅度和相位。

根据Z变换的定义

当收敛域包含|z| = 1时绝对可和，这是系统是稳定的。要保证在单位圆内收敛，需要极点都在单位圆内。(这里需要进一步理解一下收敛域相关)下面将Z变换通过单位延时转到差分方程，从时域输入输出的关系来看一下系统稳定性的解释。

即

Y(x)/X(z) = 1/(1-2(z^-1))  —>. Y(z) = 2Y(z)(z^-1) + X(z) —> y(n)=2y(n-1)+x(n)

可以看到这个系统在Z域，极点为实轴上的2，是不稳定的。从差分方程上看，每次的输入等于本次的输入加上上一次输出的2倍，这样会不断方法输入，所以会不断发散，系统不稳定。

2.5系统的频率响应

前面提到了LTI系统的特征值是复正弦信号，最终输出y(n)等于λ乘以复正弦输入信号，λ为

这里是不是和Z变换的定义有些相似？一般将频率响应写为更一般的式子：

这里就是计算系统h(n)频率响应的公式，也就是对h(n)做离散时间傅里叶变换(DTFT)。注意这里的n是从0~无穷的，也就是h(n)是一个无穷的序列，针对这种序列的傅里叶变换是离散时间傅里叶变换(DTFT)，要区别于在实际工程中更常见的离散傅里叶变换(DFT)。可见这里就是把z换成了e^jw，也就是说频率响应就是z变换H(z)在单位圆z=e^jw处的取值。由傅里叶变换得到频率响应有着明确的物理意义，但是傅里叶变换收敛条件苛刻，时域内绝对可和傅里叶变换才存在，这样导致有些信号傅里叶变换不存在。对于这种情况，拉普拉斯变换先为信号乘一个e^-σt，这样使得信号快速衰减，之后再做傅里叶变换。Z变换相当于离散域中的拉普拉斯变换，z与s的关系为z=e^s*Ts，Ts为采样周期。频率响应得到的结果一般用复数表示

也可以用幅度相位表示

|H(e^jw)|称为幅频响应，表示频率与幅度的关系。对于某些系统，在某些有些频率的信号会被抑制，表征了系统对频率的选择性。

幅频响应是带有周期性的，因为e^jw是周期旋转的，周期为2Pi。对于实系数的h(n)，幅频响应都具有偶对称的特点，即|H(e^jw)| = |H(e^-jw)|。（是不是因为只有实部，这样从0~Pi和从Pi~2Pi旋转时投射到实部的长度对称）工程中幅频响应常以db作为单位。

ψ(w)称为相频响应，也就是传递函数实部虚部间夹角。相频响应代表了输出对于相位的移动，对于不同频率的信号通过系统的时间是不一样的。另外注意的一点是输出的相位是模糊的，假如输入相位差是Pi/4，其实真实值有可能是Pi/4 + 2Pi*n的任何值。

2.5.1 群延时

直观理解为w频率附近相频响应的变化率(相位变化随频率变化的快慢)，由于是w频率附近频率，所以称为”群”。如果系统对各个频率分量有相同的延时，信号经过该系统后波形与原先一样，有一定的延时。反之如果系统对不同频率延时不同，则经过系统输出的信号将发生形变。

相频响应反应系统对输入信号延时的相对值，群延时反应系统对输入信号延时的绝对值。(如何进一步解释)如果群延时为常数，即ψ(w)=-wn，称这样的系统为线性相位。

2.6 向量的角度分析LTI系统

这里需要将传递函数由上文提到的零极点形式做变换，每项变为零极点到单位圆的向量。如下图，向量OE代表e^jw，也就是复平面上的单位圆，A、B为系统的两个零点，C、D为系统两个极点，分别可以用向量OA、OB和OC、OD表示。这样可以根据向量运算求出这些点到单位圆间的向量。

那么这些向量代表什么呢？如何与前面说的频率响应还有Z变换相关联呢？下面回顾下Z变换传递函数

2.6.1零极点形式

变换为向量形式，即零点极点到z的距离，对于频率响应而言就是零点极点到单位圆的向量的乘积。

将Z变换在单位圆上取值，代换为频率响应

可见这里每个相乘项都变成了单位圆与零点pk和极点zk的向量之差。在看本小节开头的图，对应上式，以A点举例，e^jw-zk就相当于向量AE，由于e^jw是绕着单位圆旋转，这个是以A点固定，绕单位圆旋转，用公式可以表示为

这样重新代换一下频率响应的公式，将每项代换成向量

得到幅频响应与相频响应分别如下

利用向量可以直观看出零极点位置对于频率响应的影响。幅频响应可以看作是零点到单位圆的向量积与极点到单位圆向量积的比值( (z-zk) (z-pk) 都可以看作是向量)。相频响应可以表示为 (极点个数M-零点个数N)*w+所有零点角度的和-所有极点角度的和下面看一个简单系统传函，体会一下频率响应与极点关系。

将z代换为e^jw，也就是绕单位圆转，而0.9这个点与单位圆向量的倒数，就是抚平相应。这里取0，Pi/4，Pi/2，3Pi/4，Pi等三个特殊值。

上图可以得到以下几个结论:越靠近极点，分母越小，倒数越大，幅频响应变化越快。a图上看如果逆时针旋转，超过Pi后，幅频响应对称变小，b图上来及看就是关于0对称。

2.7 两种特殊LTI系统

下面简单总结一下全通系统和最小相位系统。

2.7.1全通系统

全通系统是系统频率响应对所有频率w均为常数。最简单的全通系统为h(n)=δ(n-n0)，传递函数为H(z)=z^-n0。传递函数极点为0，这样极点到单位圆距离均为1，这样的极点对幅频响应没有影响。在实际系统中，如果想让一个系统变为全通系统，可以增加零点使之抵消掉所有的极点。实际生产中，经过一个系统可能导致相位失真，这时就需要一个幅度不变，可以调整相位的系统，将失真的相位补回来，这时就需要一个全通系统来做。假定一个系统，系数相同但排列顺序相反

可以写成

由上面两式，如果Pk为系统的极点，那么1/Pk则为系统的零点，即零点极点关于单位圆共轭倒置，这里我理解就是在复平面上互为倒数。这样互为共轭倒置的零极点就可以相互抵消。这样把z代换，系统频率响应为(这个计算是怎么来的？为什么平方？)

下面具体看一下零极点互为共轭倒置为什么可以相互抵消。这两张图可以代表全通系统零极点分布的情况，a图就是零极点互为倒数的情况，b为零极点互为共轭复数的情况。这里V1和U1并不相等。A中系统传递函数为

为了用向量分析频率响应，将上式改写为

2.7.2最小相位系统

群延时表征了系统对输入信号的延时。就像去买火车票，谁都希望去售票厅就能拿到票，而不希望排半天队，反应在相位响应上就是希望相位响应的变化越小越好。这样的系统称为最小相位系统。从零极点来看，就是所有零点都处于单位圆内的系统。为什么是零点都在单位圆内就可以呢，这里需要看一下零点对系统的影响。

对于图a的系统

用向量分析频率响应，变换为如下

根据上文从向量的观点看频率响应得出的结论，相频响应为 (极点个数M-零点个数N)*w+所有零点角度的和-所有极点角度的和。对于在单位圆内的b，w转一圈，单位圆内零点相位变化为2Pi，单位圆外零点导致相位变化为0。z^-1对系统影响为-2Pi，这样单位圆内零点导致的相位变化为0。对于更一般的情况，系统由M个极点和N个零点，对于稳定系统，所有极点都在单位圆内。n1个零点在单位圆内，n2个零点在单位圆内，