下载文档原格式
求解系统的状态⽅程求解系统的状态⽅程⼀、实验设备PC计算机，MATLAB软件，控制理论实验台⼆、实验⽬的(1)掌握状态转移矩阵的概念。
学会⽤MATLAB求解状态转移矩阵(2)学习系统齐次、⾮齐次状态⽅程求解的⽅法，计算矩阵指数，求状态响应；(3)通过编程、上机调试，掌握求解系统状态⽅程的⽅法，学会绘制输出响应和状态响应曲线；(4)掌握利⽤MATLAB导出连续状态空间模型的离散化模型的⽅法。
三、实验原理及相关基础(1)参考教材P99~101“3.8利⽤MATLAB求解系统的状态⽅程”(2)MATLAB现代控制理论仿真实验基础（3）控制理论实验台使⽤指导四、实验内容(1)求下列系统矩阵A对应的状态转移矩阵（a）（b）代码：syms lambdaA=[lambda 0 0;0 lambda 0;0 0 lambda];syms t;f=expm(A*t)（c）代码：syms t;syms lambda;A=[lambda 0 0 0;0 lambda 1 0;0 0 lambda 1;0 0 0 lambda];f=expm(A*t)(2) 已知系统a) ⽤MATLAB求状态⽅程的解析解。
选择时间向量t，绘制系统的状态响应曲线。
观察并记录这些曲线。
（1）代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];u=1;syms t;f=expm(A*t);%状态转移矩阵x0=0;s1=f*B*u;s2=int(s1,t,0,t)%状态⽅程解析解f=expm(A*t);X0=[1;0];t=[0:0.5:10];for i=1:length(t);g(i)=double(subs(f(1),t(i)));endplot(t,g)（3）状态转移矩阵syms lambdaA=[lambda 0 0;0 lambda 0;0 0 lambda];syms tf=expm(A*t)b) 计算系统在初始状态作⽤下状态响应和输出响应的数值解(⽤函数initial( )), 绘制系统的状态响应曲线和输出响应曲线。
连续系统的状态变量方程求解连续系统的状态变量方程求解通常采用数值方法，例如龙格-库塔法（Runge-Kutta）等。
在这个过程中，需要将连续系统的状态方程离散化，即将连续时间步长的微分方程转化为离散时间步长的离散方程。
求解离散方程可采用递推的方式，根据系统的初始条件和上一时刻的状态变量值，计算出当前时刻的状态变量值。
以下是一个求解连续系统状态变量方程的步骤：1. 确定连续系统的状态变量方程。
例如，给定线性定常系统dx/dt = Ax + Bu，其中x为状态变量，A和B为系统矩阵。
2. 离散化。
将状态变量方程转化为离散方程。
常见的离散化方法有前项差分变换、后项差分变换和Tustin变换。
具体变换方法取决于系统的特性以及所需的数值稳定性和精度。
例如，使用Tustin变换将连续系统离散化，得到离散状态方程x[k+1] = A*x[k] + B*u[k]。
3. 初始化。
给定初始条件，如x[0] 和u[0]，初始化状态变量值。
4. 数值求解。
使用数值方法（如龙格-库塔法）递推计算离散方程，得到一系列状态变量值x[1], x[2], ...，以及对应的输出值y[1], y[2], ...。
5. 分析结果。
根据求解得到的状态变量值和输出值，分析系统的性能，如稳定性、收敛速度等。
在MATLAB中，可以使用ode45等函数求解连续系统的状态变量方程。
以下是一个简单的示例：```MATLAB定义系统矩阵A、B和输入信号uA = [1 0; -1 1];B = [0 1];u = [1; 0.5];定义初始条件x0 = [1; 2];设置求解参数tspan = [0, 10];options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);求解状态变量方程[x, u] = ode45(@(t, x) A*x + B*u, tspan, x0, options);绘制状态变量曲线figure;plot(t, x(:, 1), 'b', 'LineWidth', 2);hold on;plot(t, x(:, 2), 'r', 'LineWidth', 2);xlabel('Time');ylabel('State Variables');legend('x1', 'x2');```这个示例中，我们使用ode45函数求解了一个线性定常系统在给定输入信号下的状态变量演化。
第三章控制系统状态方程求解3－1 线性连续定常齐次方程求解所谓齐次方程解，也就是系统的自由解，是系统在没有控制输入的情况下，由系统的初始状态引起的自由运动，其状态方程为：………………………………………………………(3－1)上式中，X是n×1维的状态向量，A是n×n的常数矩阵。
我们知道，标量定常微分方程的解为： (3)2〕与〔3－2〕式类似，我们假设〔3－1〕的解X(t)为时间t的幂级数形式，即：………………………………(3－3) 其中为与X〔t〕同维的矢量。
将〔3－3〕两边对t求导，并代入〔3－1〕式，得：上式对任意时间t都应该成立，所以变量t的各阶幂的系数都应该相等，即：即:……………………………………………〔3－4〕将系统初始条件代入〔3－3〕，可得。
代入〔3－4〕式可得：…………………………………………………………………〔3－5〕代入〔3－3〕式可得〔3－1〕式的解为： (3)6)我们记： (3)7〕其中为一矩阵指数函数，它是一个n×n的方阵。
所以〔3－6〕变为：……………………………………………………………………〔3－8〕当〔3－1〕式给定的是时刻的状态值时，不难证明：………………………………………………………………〔3－9〕从〔3－9〕可看出，形式上是一个矩阵指数函数，且也是一个各元素随时间t变化的n×n矩阵。
但本质上，它的作用是将时刻的系统状态矢量转移到t时刻的状态矢量，也就是说它起到了系统状态转移的作用，所以我们称之为状态转移矩阵(The State Transition Matrix)，并记：……………………………………………………………〔3－10〕所以：【例3－1】，求解：根据〔3－7〕式，3－2 的性质及其求法性质1：【证】根据的定义式〔3－7〕，【证毕】性质2：①②③【证】：①：根据(3－7)式，即有：②：由性质1及其关系①，③:由②式两边同时左乘，注意本身是一个n×n的方阵，，所以：即：从上式可知，矩阵指数函数的逆矩阵始终存在，且等于。
求解系统的状态方程一、实验设备PC计算机，MATLAB软件，控制理论实验台二、实验目的(1)掌握状态转移矩阵的概念。
学会用MATLAB求解状态转移矩阵(2)学习系统齐次、非齐次状态方程求解的方法，计算矩阵指数，求状态响应；(3)通过编程、上机调试，掌握求解系统状态方程的方法，学会绘制输出响应和状态响应曲线；(4)掌握利用MATLAB导出连续状态空间模型的离散化模型的方法。
三、实验原理及相关基础(1)参考教材P99~101“3.8利用MATLAB求解系统的状态方程”(2)MATLAB现代控制理论仿真实验基础（3）控制理论实验台使用指导四、实验内容(1)求下列系统矩阵A对应的状态转移矩阵（a）（b）代码：syms lambdaA=[lambda 0 0;0 lambda 0;0 0 lambda];syms t;f=expm(A*t)（c）代码：syms t;syms lambda;A=[lambda 0 0 0;0 lambda 1 0;0 0 lambda 1;0 0 0 lambda];f=expm(A*t)(2) 已知系统a) 用MATLAB求状态方程的解析解。
选择时间向量t，绘制系统的状态响应曲线。
观察并记录这些曲线。
（1）代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];u=1;syms t;f=expm(A*t);%状态转移矩阵x0=0;s1=f*B*u;s2=int(s1,t,0,t)%状态方程解析解状态曲线：（2）A=[0 1;-2 -3];syms t;f=expm(A*t);X0=[1;0];t=[0:0.5:10];for i=1:length(t);g(i)=double(subs(f(1),t(i)));endplot(t,g)（3）状态转移矩阵syms lambdaA=[lambda 0 0;0 lambda 0;0 0 lambda];syms tf=expm(A*t)b) 计算系统在初始状态作用下状态响应和输出响应的数值解(用函数initial( )), 绘制系统的状态响应曲线和输出响应曲线。
观察并记录这些响应曲线，然后将这一状态响应曲线与a)中状态响应曲线进行比较。
代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];G=ss(A,B,C,D);t=[0:0.5:10];x0=[1;0][y0,t,x0]=initial(G,x0,t); plot(t,x0,'-',t,y0,'-')c) 根据b)中所得的状态响应的数值解，绘制系统的状态轨迹(用命令plot(x(:,1), x(:,2)))。
记录系统状态转移的过程，结合a)和b)中的状态响应曲线分析这一过程。
代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];t=[0:0.01:10];x0=[1;0];G=ss(A,B,C,D)[y,t,x]=initial(G,x0,t);plot(x(:,1),x(:,2))2) 令初始状态为零，输入为u(t)=1(t).a)用MATLAB求状态方程的解析解。
选择时间向量t，绘制系统的状态响应曲线。
观察并记录这些曲线。
代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];G=ss(A,B,C,D);[y,t,x]=step(G);plot(t,x)b)计算系统在初始状态作用下状态响应和输出响应的数值解, 绘制系统的状态响应曲线和输出响应曲线。
观察并记录这些响应曲线，然后将这一状态响应曲线与a).中状态响应曲线进行比较。
代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];G=ss(A,B,C,D);G=ss(A,B,C,D);t=[0:0.5:10];x0=[1;-1];[y0,t,x0]=initial(G,x0,t);plot(t,x0,'-',t,y0,'-')c) 根据b)中所得的状态响应的数值解，绘制系统的状态轨迹。
记录系统状态转移的过程，结合a)和b)中的状态响应曲线分析这一过程。
代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];t=[0:0.5:10];G=ss(A,B,C,D);x0=[0 0];[y0,t,x0]=initial(G,x0,t);plot(t,x0,'-',t,y0,'-')绘制系统的状态响应曲线、输出响应曲线和状态轨迹。
观察和分析这些响应曲线和状态轨迹是否是（1）和（2）中的响应曲线和状态轨迹的叠加。
代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];t=[0:0.01:10];x0=[1 -1];G=ss(A,B,C,D);[y,t,x]=initial(G,x0,t);plot(t,x)4）令初始状态为零，输入为u(t)=3sin(5t)。
计算状态响应和输出响应的数值解(用函数lsim( ))，并绘制系统的状态响应曲线、输出响应曲线和状态轨迹。
代码：A=[0 1; -2 -3];B=[3;0];C=[1 1];D=[0];t=[0:0.01:10];u=3*sin(5*t);G=ss(A,B,C,D);[y,t,x]=lsim(G,u,t); plot(t,x)（3）已知系统1）当输入为u(t)=(t)时，用函数initial( )和 impulse( )求解系统的状态响应和输出响应的数值解，并绘制系统的状态响应曲线、输出响应曲线和状态轨迹。
状态响应：A=[0,1,0;0,0,1;-6,-11,-6];B=[0;0;1];C=[6,0,0];D=0;t=[0:0.5:10];G=ss(A,B,C,D);x0=[1,0,-1];[y,t,x]=initial(G,x0,t);u=ones(size(t));plot(t,x,t,y)输出响应：A=[0,1,0;0,0,1;-6,-11,-6]; B=[0;0;1];C=[6,0,0];D=0;t=[0:0.01:10];u=ones(size(t));G=ss(A,B,C,D);[y,t,x]=lsim(G,u,t);plot(t,x)2）当输入为u(t)=1(t)时，用函数initial( )和 step( )求解系统的状态响应和输出响应的数值解，并绘制系统的状态响应曲线、输出响应曲线和状态轨迹。
状态响应A=[0,1,0;0,0,1;-6,-11,-6];B=[0;0;1];C=[6,0,0];D=0;t=[0:0.5:10];G=ss(A,B,C,D);x0=[1,0,-1];[y,t,x]=initial(G,x0,t);u=step(G);plot(t,x,t,y)输出响应：A=[0,1,0;0,0,1;-6,-11,-6]; B=[0;0;1];C=[6,0,0];D=0;G=ss(A,B,C,D);[y,t,x]=step(G);plot(t,x)3）当输入为u(t)= t 时，用函数initial( )和 lsim( )求解系统的状态响应和输出响应的数值解，并绘制系统的状态响应曲线、输出响应曲线和状态轨迹。
状态响应A=[0,1,0;0,0,1;-6,-11,-6];B=[0;0;1];C=[6,0,0];D=0;t=[0:0.5:10];G=ss(A,B,C,D);x0=[1,0,-1];[y,t,x]=initial(G,x0,t);u=t;plot(t,x,t,y)输出响应A=[0,1,0;0,0,1;-6,-11,-6]; B=[0;0;1];C=[6,0,0];D=0;t=[0:0.01:10];u=t;G=ss(A,B,C,D);[y,t,x]=lsim(G,u,t);plot(t,x)4）当输入为时，用函数initial( ) 和lsim( )求解系统的状态响应和输出响应的数值解，并绘制系统的状态响应曲线、输出响应曲线和状态轨迹.状态响应A=[0,1,0;0,0,1;-6,-11,-6];B=[0;0;1];C=[6,0,0];D=0;t=[0:0.5:10];G=ss(A,B,C,D);x0=[1,0,-1];[y,t,x]=initial(G,x0,t);u=sin(t);plot(t,x,t,y)输出响应A=[0,1,0;0,0,1;-6,-11,-6]; B=[0;0;1];C=[6,0,0];D=0;t=[0:0.01:10];u=sin(t);G=ss(A,B,C,D);[y,t,x]=lsim(G,u,t);plot(t,x)(4) 已知一个连续系统的状态方程是若取采样周期秒0.05 T1）试求相应的离散化状态空间模型；代码:syms T;A=[0 1;-25 -4];B=[0;1];[Gz,Hz]=c2d(A,B,T)G =[ exp(- 2*T - 21^(1/2)*T*i)/2 + exp(- 2*T + 21^(1/2)*T*i)/2 + (21^(1/2)*exp(-2*T - 21^(1/2)*T*i)*i)/21 - (21^(1/2)*exp(- 2*T + 21^(1/2)*T*i)*i)/21, (21^(1/2)*exp(- 2*T - 21^(1/2)*T*i)*i)/42 - (21^(1/2)*exp(- 2*T +21^(1/2)*T*i)*i)/42][ - (21^(1/2)*exp(- 2*T -21^(1/2)*T*i)*25*i)/42 + (21^(1/2)*exp(- 2*T + 21^(1/2)*T*i)*25*i)/42, exp(-2*T - 21^(1/2)*T*i)/2 + exp(- 2*T + 21^(1/2)*T*i)/2 - (21^(1/2)*exp(- 2*T -21^(1/2)*T*i)*i)/21 + (21^(1/2)*exp(- 2*T + 21^(1/2)*T*i)*i)/21]H =1/25 - exp(- 2*T + 21^(1/2)*T*i)/50 - (21^(1/2)*exp(- 2*T -21^(1/2)*T*i)*i)/525 + (21^(1/2)*exp(- 2*T + 21^(1/2)*T*i)*i)/525 - exp(- 2*T- 21^(1/2)*T*i)/50(21^(1/2)*exp(- 2*T - 21^(1/2)*T*i)*i)/42 - (21^(1/2)*exp(- 2*T +21^(1/2)*T*i)*i)/422）分析不同采样周期下，离散化状态空间模型的结果。
