状态方程simulink仿真_控制系统设计与仿真作业与复习资料

仿真的三次作业和三次上机实验报告，百度云链接: http://pan.baidu.com/s/1eRsWZ2i 密码: 9s8v注意不要直接抄袭喔！如果老师发现所有学生的作业都一样的话你们懂的

当时考试的时候总结的，考题基本【或者肯定】就从这里出……

---

采用四阶龙格库塔法求如下二阶系统的阶跃响应。

解：

T=10;%仿真时长b=0.5;w=5; %%%%%%%%%%%采用矩阵方法的四阶龙格库塔法%%%%%dt=0.05;%积分步长N=T/dt;%积分总步数Y1=[0;0];YY=Y1;A=[0 1;-w^2 -2*b*w];%%%%状态观测矩阵%%%%B=[0;1];C=[1 0];U=w^2;for i=1:N;    K1=A*YY+B*U;    K2=A*(YY+dt/2*K1)+B*U;    K3=A*(YY+dt/2*K2)+B*U;    K4=A*(YY+dt*K3)+B*U;    YY=YY+dt/6*(K1+2*K2+2*K3+K4);    Y1=[Y1 YY];end;figureplot([0:dt:T],C*Y1);grid on

求脉冲响应：

%把U定义在for中for i=1:N;if i < 1/dt        U = w^2;    else        U = 0;    end    K1=A*YY+B*U;    K2=A*(YY+dt/2*K1)+B*U;    K3=A*(YY+dt/2*K2)+B*U;    K4=A*(YY+dt*K3)+B*U;    YY=YY+dt/6*(K1+2*K2+2*K3+K4);    Y1=[Y1 YY];end;

ode45解：

function r=second_order(t,x);b=0.5;w=5;r=[x(2);-2*b*w*x(2)-w^2*x(1)+w^2];endt=[0 T];x0=[0 0];[t y]=ode45('second_order',t,x0);figure;hold onplot(t,y(:,1));

---

采用四阶龙格库塔法求高阶系统单位阶跃响应曲线的数值解，并从数值解中得到系统的上升时间、峰值时间、调节时间等时域指标。

解：

t=40;%仿真时长b=0.5;w=10;T = 5; %%%%%%%%%%%采用矩阵方法的四阶龙格库塔法%%%%%dt=0.05;%积分步长N=t/dt;%积分总步数Y1=[0;0;0];YY=Y1;A=[0 1 0;0 0 1;-w^2/T (-2*b*w-T*w^2)/T -(1+2*b*w)/T];B=[0;0;1];C=[1 0 0];U = w^2/T;for i=1:N;        K1=A*YY+B*U;    K2=A*(YY+dt/2*K1)+B*U;    K3=A*(YY+dt/2*K2)+B*U;    K4=A*(YY+dt*K3)+B*U;    YY=YY+dt/6*(K1+2*K2+2*K3+K4);    Y1=[Y1 YY];end;figureplot(0:dt:t,C*Y1);

响应曲线参数

%y(1,:)是矩阵算出的(C*Y1)，y(:,1)是ode45算出的tr = t(find(y(:,1) > 1,1));%上升时间ts = t(max(find(y(:,1)<0.98|y>(:,1)1.02,1,'last')));%调节时间tp = t(y == max(y(:,1)));%峰值时间overshoot = max(y(:,1)) - 1;%超调量

---

单位反馈系统的开环传递函数G(s)=10/(s(s+1))，求闭环传递函数的多项式展开形式，零极点形式，状态空间形式以及部分分式形式。

解：

clear all;num = 10;den = [1 1 0];GH = tf(num,den);%G(s)H(s)G1= feedback(GH,1)%G1为单位闭环负反馈[numc,denc] = cloop(num,den);[z,p,k] = tf2zp(numc,denc)%将系统转化成零极点形式[A,B,C,D] = tf2ss(numc,denc)%将系统转化成状态方程形式[r,p,k] = residue(numc,denc)%将系统转化成部分分式展开形式

---

大意就是把下面的非线性方程组在Simulink里用框图的格式画出来。

参数入口为α，β，γ的值以及xi，的初值。(注：α = 8/3，β = 10，γ = 28，初值分别为x1 = 0，x2 = 0，x3 = 0.001时，系统输出为混沌状态)

解：

---

判断如下系统的稳定性：

解：

num = [1 2 3 1];den = [1 5 2 1 1];%第一种方法：检查特征多项式的根r = roots(den);i1 = find(real(r) > 0);n1 = length(i1);if(n1 > 0)    disp('The System is Unstable!');else    disp('The System is Stable!');end%第二种方法：检查极点[z,p,k] = tf2zp(num,den);i2 = find(real(p) > 0);n2 = length(i2);if(n2 > 0)    disp('The System is Unstable!');else    disp('The System is Stable!');end%第三种方法：零极点图pzmap(num,den);title('Zero-Pole Map');%第四种方法：频域margin(num,den);

解：

clear all;num = [-3 2];den = [1 -0.2 -0.25 0.05];%第一种方法：检查特征多项式的根r = roots(den);i1 = find(abs(r) > 1);n1 = length(i1);if(n1 > 0)    disp('The System is Unstable!');else    disp('The System is Stable!');end%第二种方法：检查极点[z,p,k] = tf2zp(num,den);i2 = find(abs(r) > 1);n2 = length(i2);if(n2 > 0)    disp('The System is Unstable!');else    disp('The System is Stable!');end%第三种方法：零极点图pzmap(num,den);title('Zero-Pole Map');

解：

clear all;A = [1 3;5 2];%第一种方法：检查矩阵特征值r1 = eig(A);i1 = find(real(r1) > 0);n1 = length(i1);if(n1 > 0)    disp('The System is Unstable!');else    disp('The System is Stable!');end%第二种方法：检查矩阵的特征多项式p = poly(A);r2 = roots(p);i2 = find(real(r2) > 0);n2 = length(i2);if(n2 > 0)    disp('The System is Unstable!');else    disp('The System is Stable!');end%第三种方法：李雅普诺夫第二法Q = eye(size(A));P = lyap(A,Q);i21 = find(P(1,1) > 0);n21 = length(i21);i22 = find(det(P) > 0);n22 = length(i22);if(n21 > 0 && n22 > 0)    disp('The System is Stable!');else    disp('The System is Unstable!');end

---

试产生一周期为5秒，时长为30秒，最大值为1，最小值为0的三角波，并分析如下一阶系统在三角波输入下的时间响应曲线。

G(s)=1/(s+1)

解：

clear all;t = 0:0.001:30;k = length(t);y = zeros(1,k);for i = 1:k    if (t(i) >= 0 && t(i) <= 2.5)        y(i) = t(i)/2.5;    elseif (t(i) > 2.5 && t(i) <= 5)        y(i) = 2 - t(i)/2.5;     else        y(i) = y(i - (k-1)/6);    endend        num = 1;den = [1 1];lsim(num,den,y,t)

---

已知系统开环传递函数：

1)试用根轨迹方法求解其临界稳定增益k。

2)若k = 10，试绘制伯德图，求相角裕度和幅值裕度，并判断其稳定性。

解：

clear all;num = 1;den = conv([2 1],conv([1 1],[0.1 1]));rlocus(num,den);figure(2);margin(10,den);

---

已知系统开环传递函数：

试设计超前校正环节，使校正后系统的静态速度误差系数Kv≤6，相角裕度为45°，并绘制校正前后系统的单位阶跃响应曲线和开环Bode图加以比较。

解：

clear all;%(1)求K0%因为Kv<=6，所以K0 = 3;%(2)作原系统的波特图与阶跃响应曲线，检查是否满足题目要求n1 = 2;d1 = conv(conv([1 0],[0.1 1]),[0.3 1]);[mag1,phase1,w1] = bode(K0*n1,d1);figure(1);margin(mag1,phase1,w1);hold on;figure(2);s1 = tf(K0*n1,d1);sys1 = feedback(s1,1);step(sys1);hold on;%(3)求超前校正器的传递函数[mag,phase,w] = bode(s1);gama = 55;gam = gama*pi/180;alfa = (1-sin(gam))/(1+sin(gam));adb = 20*log10(mag);am = -10*log10(1/alfa);wc = spline(adb,w,am);T = 1/(wc*sqrt(alfa));Gc = tf([T 1],[alfa*T 1]);%(4)校验系统校正后是否满足要求，并给出系统校正后的阶跃响应曲线s2 = s1*Gc;[mag2,phase2,w2] = bode(s2);figure(1);margin(mag2,phase2,w2);figure(2);sys2 = feedback(s2,1);step(sys2);

---

已知系统开环传递函数：

试设计超前校正环节，使得校正后系统的闭环主导极点满足：单位阶跃响应超调量为20%，调节时间为1s，单位斜坡响应的稳态误差为10%。

解：

clear all;%(1)求校正器增益%因稳态误差为0.1，则2Kc=10，Kc=5%(2)校验原系统是否满足要求t=0:0.01:20;d1 = conv(conv([1 0],[0.1 1]),[0.3 1]);sa = tf(2,d1);sysa = feedback(sa,1);[y,t] = step(sysa,t);overshoot = max(y) - 1;ts = t(max(find(y<0.98|y>1.02,1,'last')));%(3)确定期望极点的位置overshoot0 = 0.2;s = ((log(1/overshoot0))^2/((pi)^2+(log(1/overshoot0))^2))^(1/2);wn = 4.5/s;p = [1 2*s*wn wn*wn];rt = roots(p);%(4)求校正器的传递函数kc=5;s_1=rt(1);nk1=2;dk1=conv(conv([1 0],[0.1 1]),[0.3 1]);ngv=polyval(nk1,s_1);dgv=polyval(dk1,s_1);g=ngv/dgv;zetag=angle(g);mg=abs(g);ms=abs(s_1);zetas=angle(s_1);tz=(sin(zetas)-kc*mg*sin(zetag-zetas))/(kc*mg*ms*sin(zetag));tp=-(kc*mg*sin(zetas)+sin(zetag+zetas))/(ms*sin(zetag));nk=[tz,1];dk=[tp,1];Gc=tf(nk,dk);%(5)校验校正器n1 = 10;s1 = tf(n1,d1);sys1 = feedback(s1*Gc,1);[y1,t] = step(sys1,t);overshoot1 = max(y1) - 1;ts1 = t(max(find(y1<0.98|y1>1.02,1,'last')));%(6)校正前后对比step(sysa);hold on;step(sys1);

---

已知单位负反馈系统的开环传递函数：

对系统采用比例—积分(PI)控制，比例系数为Kp = 2，试求积分时间常数为TI = 3,6,14,21,25时，闭环系统的单位阶跃响应曲线，分析时间常数对控制效果的影响。

解：

clear all;G = tf(1,conv([1 1],[2 1]));Kp = 2;Ti = [3 6 14 21 25];for i = 1:5    G1 = tf([Kp Kp/Ti(i)],[1 0]);    sys = feedback(G1*G,1);    step(sys);    hold on;end