实验六 离散系统状态方程的求解

实验六 离散系统状态方程的求解

一、实验目的

（1）了解离散系统状态方程求解方法。

（2）了解离散系统信号流图化简的方法。

（3）了解函数ode45的调用方法。

二、实验原理

离散系统状态方程的一般形式为 x(k+1)=Ax(k)+Bf(k)

在些只对单输入的n 阶离散系统的状态议程求解。一般采用递推迭代的方式求解，由裙的条件x(0)和激励f(0)求出k=1时的x(1)，然后依次迭代求得所要求的x(0)，……，x(n)的值。

编程时应注意,MATLAB 中变量下标不允许为零，则裙点的下标只能取1，第n 步的x 的下标为n+1。

三、涉及的MATLAB 函数

zeros(2,1)

y=lsim(sys,f,[],x0)

for i=1:n end

clear all

采用函数ode45可以求解微分方程。其调用格式如下

[t,y]=ode45(odefun,tspan,y0)

其中，odefun 指状态方程的表达式，tspan 指状态方程对应的起止时间]t0,tf]，y0指状态变量的初始状态。

四、实验内容与方法

1．验证性实验

采用MATLAB 语言编程，求解离散系统状态方程，并绘制状态变量的波形。

（1）已知离散系统的状态方程为 )(01)1()1(25.025.005.0)1()1(2121k f k x k x k x k x ⎥⎦

⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡++⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡++ 初始条件为x(0)=⎥⎦

⎤⎢⎣⎡-5.01,激励为f(k)=0.5ε(k)，确定该状态方程x(k)前10步的解，并画出波形。

MATLAB 程序：

%离散系统状态求解

%A=input(‘系统矩阵 A=’)

%B=input(‘系数矩阵 B=’)

%x0=input(‘初始状态矩阵 x0=’)

%n=input(‘要求计算的步长n=’)

%f=input(‘输入信号f=’) %需求长度为n的数组

clear all

A=[0.5 0;0.25 0.25];

B=[1;0];

x0=[-1;0.5];

n=10;

f=[0 0.5*ones(1,n-1)];

x(:,1)=x0;

for i=1:n

x(:,i+1)=A*x(:,i)+B*f(i);

end

subplot(2,1,1);stem([0:n],x(1,:));

subplot(2,1,2);stem([0:n],x(2,:));

离散系统状态方程的求解结果如图6.1所示

图6－1 离散系统状态方程的求解（2）离散系统状态求解

MATLAB程序：

A=[0 1;-2 3];B=[0;1];

C=[1 1;2 -1];D=zeros(2,1);

x0=[1;-1];%初始条件

N=10;f=ones(1,N);

sys=ss(A,B,C,D,[]);

y=lsim(sys,f,[],x0);k=0:N-1;

subplot(2,1,1);

stem(k,y(:,1),'b');

subplot(2,1,2);stem(k,y(:,2),'b');

离散系统状态议程的求解结果如图6.2所示。

图6－2 离散系统状态方程的求解

2.程序设计实验

(1)离散系统状态方程为： x(k+1)=Ax(k)+Bf(k)

其中A=⎥⎦⎤⎢⎣⎡25.025.005.0，B=⎥⎦

⎤⎢⎣⎡01，初始状态⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡00)0()0(21x x ，激励f(k)=δ(k)，确定该状态方程x(k)前10步的解，并画出波形。

MATLAB 程序：

%离散系统状态求解

%A=input(‘系统矩阵 A=’)

%B=input(‘系数矩阵 B=’)

%x0=input(‘初始状态矩阵 x0=’)

%n=input(‘要求计算的步长 n=’)

%f=input(‘输入信号 f=’) %需求长度为n 的数组

clear all

A=[0.5 0;0.25 0.25];

B=[1;0];

x0=[0;0];

n=10;

f=[1 zeros(1,n-1)];

x(:,1)=x0;

for i=1:n

x(:,i+1)=A*x(:,i)+B*f(i);

end

subplot(2,1,1);stem([0:n],x(1,:));

subplot(2,1,2);stem([0:n],x(2,:));

求解结果如图6.3所示：

图6－3 离散系统状态议程的求解

（2）描述离散时间系统的信号流图如图所示，确定该系统的系统函数（离散系统信号流图的形式与连续系统相同，只不过是变量s换为z，在此不再详述。）。

MATLAB程序：

syms z;

Q=[0 0 0 0 0 0;2 0 -3 0 -2 0;0 1/z 0 0 0 0;3 0 1 0 0 0;0 0 0 1/z 0 0;1 0 0 0 2 0];

B=[1;0;0;0;0;0];

I=eye(size(Q));

H=(I-Q)\B;

H6=H(6);

pretty(H6);

执行后得到结果：

2

24 + 9 z + z

-------------

2

z + 3 z + 2