机械工程控制理论课程作业一、计算题(应用MATLAB 求解)1. 一系统由下列两个子系统并联而成，试确定该系统的整体状态方程模型、传递函数模型，并确定系统的零、极点。
如取采样周期T=0.1s ，确定该系统所对应的Z 传递函数和离散状态方程，并判别系统的稳定性。
子系统1：系统状态空间模型的参数[]0,21,01,0152==⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡--=D C B A子系统2：系统的传递函数模型为33)(+=s s G答： Matlab 程序： clc;clear all;A=[-2,-5;1,0];B=[1;0];C=[1,2];D=0;T=0.1;ss1=ss(A,B,C,D); %子系统1的状态空间方程[num1,den1]=ss2tf(A,B,C,D); %子系统1传递双数的分子分母各阶系数 sys1=tf(num1,den1); %子系统1转化为传递函数模型 sys2=tf([3],[1,3]); %子系统2的传递函数模型ssa=ss1+sys2; %两系统并联，系统总体状态空间方程 sys=sys1+sys2; %两系统并联，系统总体传递函数 zero(sys) %系统的零点 pole(sys) %系统的极点 ssad=c2d(ssa,T) %系统离散状态方程 sysad=c2d(sys,T) %系统脉冲传递函数if max(abs(pole(sysad)))>1 %判断是否稳定，通过极点我们可以判断系统稳定性 disp('系统不稳定') elsedisp('系统稳定') end2. 时不变系统 Cx y Bu Ax x =+=.，且 ⎥⎦⎤⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=214321,001001,716531313C B A , 试计算该系统的特征值，并判别其能控性与能观性，确定系统状态方程模型(对角标准型)Jordan 标准型。
答： Matlab 程序： clc;clear all;A=[-3,1,3;1,-3,5;6,1,7]; B=[1,0;0,1;0,0;]; C=[1,2,3;4,1,2];D=0;sys=ss(A,B,C,D) %建立状态方程模型 e = eig(sys) ; %求系统特征值 rank(ctrb(A,B)) %系统能控性判断 rank(obsv(A,C)) %系统能观性判断[csys,T]=canon(sys,'modal')%化为对角标准型,T 为转化矩阵3. 若系统的状态方程模型参数⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=5.34231201,0012321301000010B A 选择加权矩阵Q=diag{1, 2, 3, 4}(diag = diagonal matrix 对角矩阵)及R=eye(2)(identity matrix 单位矩阵), 则设计出这一线性二次型指标的最优控制器及在最优控制下的闭环系统极点位置。
答： Matlab 程序： clc;clear all;A=[0 1 0 0;0 0 1 0;-3 1 2 3;2 1 0 0];B=[1 0;2 1;3 2;4 3.5];C=eye(4);D=zeros(4,2); sys1=ss(A,B,C,D);disp('开环系统极点位置')eig(A) %开环系统极点位置 figure(1);step(sys1); %查看系统状态的响应 Q=diag([1:4]);R=eye(2);[K,S,E] = LQR(A,B,Q,R); %线性二次型指标的最优控制器 sys2=ss(A-B*K,B,C,D); %加控制器后的反馈系统 disp('闭环系统极点位置')eig(A-B*K) %闭环系统极点位置 figure(2);step(sys2); %加控制器后系统阶跃响应4. 已知线性离散系统的状态方程，试判断系统的稳定性： (a) )(05.05.01)(kT x T kT x ⎥⎦⎤⎢⎣⎡=+, (b) )(368.0632.0)(632.0632.0632.0368.0)(kT u kT x T kT x ⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡-=+答:Matlab 程序： clc;clear all;A1=[1 0.5;0.5 0];A2=[0.368 0.632;-0.632 0.632]; disp('A1特征根的模')e1=eig(A1); %求A1系统的特征值f1=find(abs(e1)>1) %所以特征值中模大于1的个数 if f1>0 %判断系统稳定性 disp('A1不系统稳定') elsedisp('A1系统稳定') enddisp('A2特征根的模')e2=eig(A2); %同理 f2=find(abs(e2)>1) if f2>0disp('A2系统不稳定') elsedisp('A2系统稳定') end二、编程题(选做2题)1. 现有一组开环系统频率特性数据（G (jw i )H (jw i ), i =1,2…,n ），试编写一小程序，来判断对应的闭环系统是否稳定.function exe21 B=[1];A=[1 2 1 2 1]; %假设试验传递函数[H,W]=freqs(B,A); %试验产生开环的频率相应数据%%%%由试验数据估计出系统开环的传递函数 disp('试验数据GH(Wj)：') Hsum0=sum(H); for b=0:50for a=(b+1):51[B,A]=invfreqs(H,W,b,a); H1=freqs(B,A,W);sum1=sum(H1);if (norm((sum1)-(sum0)))<1e-5 %比较试验数据和估计数据之间的误差G=tf(B,A); break endendif (norm((sum1)-(sum0)))<1e-5 %比较试验数据和估计数据之间的误差G=tf(B,A); break endendG; %开环传递函数GG=1+G; %闭环传递函数G/(1+G*H)的分母[num,den]=tfdata(GG);num=num{:}; %1+G*H 的分子即为闭环传递函数的特征式a=roots(num); %1+G*H 的分子即为闭环传递函数的特征式，由此求根判断flag=0; for i=1:length(a) if real(a(i))>0 flag=flag+1; end end if flag>0disp('系统闭环不稳定') elsedisp('系统闭环稳定')end2. 现通过实验取得一线性系统的频率特性参量，相关数据存在文件FredataA.mat ，应用invfreqs 语句估计该系统参数(传递函数分子、分母系数)。
3. 编写求解线性时变系统的时域响应的子程序。
答：运用Matlab 计算线性时变系统的时域响应。
程序定义了名为time_respones3的函数，函数的形参有系数矩阵A 、B 、C 、D 和输入向量u ，初始状态向量xt0、初始时刻、以及绘制曲线图的终端时间tmax 和状态转移矩阵的积分项数N 等。
状态转移矩阵通过定义的子函数matrix_exp 获得。
原理：线性时变连续系统的状态转移矩阵计算原理：（1） 幂级数法:当且仅当A (t )A (τ)=A (τ)A (t )时，状态转移矩阵为∅(t,t 0)=e ∫A(τ)dτt t 0=I +∫A(τ)dτt t 0+12!(∫A(τ)dτt t 0)2+⋯+1n!(∫A (τ)dτt t 0)n +⋯（2） 递推级数法:不需要额外的条件，∅(t,t 0)=I +∫A(τ1)dτ1tt 0+∫∫A(τ1)A(τ2)dτ2dτ1τ1t 0tt 0+∫∫∫A(τ1)A(τ2)A(τ3)dτ3dτ2dτ1τ2t 0τ1t 0tt 0+⋯系统状态受控运动原理:x (t )=∅(t,t 0)x (t 0)+∫∅(t,t 0)B(τ)u(τ)dτtt 0系统输出响应：求出状态向量x(t)后，将其代入输出方程，即可确定输出响应，即y (t )=C (t )x (t )+D(t)u(t)主程序：function [phit0,x,y] =time_response3(A,B,C,D,u,xt0,t0,tmax,N) %A,B,C,D 为系数矩阵；u 为输入量，xt0为初始状态向量 %N 为状态转移矩阵积分项数，tmax 为曲线时间上限值。
syms t tao; %定义变量 phit0=matrix_exp(A,N);phit=simplify(subs(phit0,'t0',t0));xx0=phit*xt0; %自由分量phitao=subs(phit0,'t0','tao'); %phit （t ，t0）改为phit （t ，tao ） Btao=subs(B,'t','tao'); %B(t)改为B （tao ）utao=subs(u,'t','tao'); %u(t)改为u （tao ） xu=int(phitao*Btao*utao,tao,t0,t); %受控分量 x=simplify(xx0+xu); %状态响应 y=simplify(C*x+D*u); %输出响应 q=size(C,1); for k=1:qsubplot(2,2,k);ezplot(t,y(k),[t0,tmax]);grid;title(' ') xlabel('时间 t'); ylabel('输出量 y') end end 子程序：function [phit] = matrix_exp(A,N) %A 系统矩阵；N 为积分项数syms t t0 tao; %变量I=eye(size(A));int_AA=I;phit=I;Atao=subs(A,'t','tao'); %将A 矩阵中的符号变量t 替换为新的值tao if(A*Atao==Atao*A) %应用幂级数法计算int_A=int(Atao,tao,'t0','t'); %定义积分函数Atao ，积分变量tao ，积分区间t0，t for i=1:Nint_AA=simplify(int_A*int_AA/i); %对表达式进行简化 phit=phit+int_AA;endelse %用递推法计算 for i=1:NAA=Atao*subs(int_AA,'t','tao'); if(AA==0) break; endint_AA=simplify(int(AA,tao,t0,t)); phit=phit+int_AA; end end end程序运行例子： 假设系统的动态方程为：2121002122x y u e x x e te x x t t t =⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡--- 若初始状态为x （0）=0试求u=3（t>=0）作用下的系统响应。