中国矿业大学控制系统计算机仿真实验答案

实验一 1、（教材P37页 例2.1）已知一阶微分方程：

?

y(0)?1/3,h?0.05试编写程序，用欧拉法求 y30的值；并画出图形。 y(t?1.5)，即

解：程序如下： clc clear

Tf=input('仿真时间Tf='); %输入仿真时间 h=input('计算步长h='); %输入仿真步长 x0=1/3; t=0; x=x0; %输入初始值 for i=1:Tf/h K1=-30*x0; x1=x0+h*K1;

x=[x;x1]; %变量以向量形式保存

t=[t;t(i)+h]; %对应时刻以向量形式保存 x0=x1; end

[t,x] %以数据形式输出 plot(t,x) %以曲线形式输出 所绘图形如下：

0.5y??30y,0.40.30.20.1X: 1.5Y: 3.104e-0100-0.1-0.200.20.40.60.811.21.41.6?1?2、2. （教材P79习题2.5）已知系统的状 ?x态方程和输出方程为： ?x??式中 u(t)=1(t)。初始条件为： ?2?x1(0)=x2(0)=x3(0)=0。取h=0.05， ??3?x??试用RK4法求t=0.5时的y 值。

解：程序如下：

1 0??x1??0???8 ? ?x???4?u???19 0 1???2??????12 0 0????10???x3????x1??y?[1 0 0]?x?2?1 ???x3?

clc clear u=1;

x0=0; t=0;y=0; x0=[0;0;0];

A=[-8 1 0;-19 0 1;-12 0 0]; B=[0;4;10]; C=[1 0 0];

Tf=input('仿真时间Tf='); h=input('计算步长h='); x=x0;

for i=1:Tf/h K1=A*x+B*u;

K2=A*(x+h*K1/2)+B*u; K3=A*(x+h*K2/2)+B*u; K4=A*(x+h*K3)+B*u;

x=x+h*(K1+2*K2+2*K3+K4)/6; y=[y;C*x]; t=[t;t(i)+h]; end y

plot(t,y)

所绘图形如下：

0.90.80.70.60.50.40.30.20.1005101520253035

0.250.2X: 0.5Y: 0.22010.150.10.05000.050.10.150.20.250.30.350.40.450.52

3、已知系统的传递函数， 试求其状态方程及输出方程。 （教材P26页 习题1.7 ） s?3G(s)?32s?3s?3s?1

解：程序如下： den=[1 3];

num=[1 3 3 1];

[A,B,C,D]=tf2ss(den,num) 运行结果如下： A =

-3 -3 -1 1 0 0 0 1 0 B =

1 0 0 C =

0 1 3 D =

0

4、已知线性定常系统的状态方程

??及输出方程，试求系统的传递函数。?x1???1 1 ???（教材P26页 习题1.8 ） ? x?2 0 ?1 ?x????解：程序如下：

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

B=[0;0;1];

?x1?C=[1 0 1]; ??[den,num]=ss2tf(A,B,C,D) y?[1 0 1]运行结果如下：

?x2???x3??

0???x1??0?0-2? ?x???0?u??2?????x3?????1??3

den =

0 1 2 1

num =

1 4 5 2

5、（教材P48页 例2.4 ）如图所示的直流电机拖动系统，试在SIMULINK环境下

建立系统仿真模型：

1）采用欧拉法与四阶龙格库塔法进行仿真，并确定其步长的稳定区间，取不同步长比较仿真结果；

2）采用RK4研究外环PI控制器参数的改变(0.17 0.5 1 1.5)，观其阶跃响应对系统动态参数的影响。

（1）、欧拉法 clc clear

Tf=10; %输入仿真时间

h=input('计算步长h='); %输入仿真步长 P=[1 0.01 0.1 0; 0 0.85 1 0.17; 1 0.01 1 0;

0 0.051 1 0.15; 1 0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0; 1 0.01 0.1 0;

1 0.01 0.0044 0]; %输入各环节参数 W0=zeros(9,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1; W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1; W(8,6)=1;W(9,7)=1; %输入系统连接矩阵

4

A=diag(P(:, 1));B=diag(P(:, 2));

C=diag(P(:, 3));D=diag(P(:, 4)); %生成系数矩阵 Q=B-D*W;Q1=inv(Q); R=C*W-A; V=C*W0;

A1=Q1*R; B1=Q1*V; %生成闭环系数矩阵 C1=[0 0 0 0 0 0 1 0 0];

u=1; %阶跃输入幅值

x0=zeros(9,1); t=0;y=0; %设置初值 x=x0;

for i=1:Tf/h

K1=A1*x+B1*u; x=x+h*K1; y=[y;C1*x]; t=[t;t(i)+h]; end

plot(t,y)

计算步长h=0.011 计算步长h=0.012

403530250.52.521.51x 105020015-0.51050-1-1.5-2012345678910012345678910

四阶龙格库塔法： clc clear

Tf=input('仿真时间Tf='); %输入仿真时间 h=input('计算步长h='); %输入仿真步长 P=[1 0.01 0.1 0; 0 0.85 1 0.17; 1 0.01 1 0;

0 0.051 1 0.15; 1 0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0;

5

1 0.01 0.1 0;

1 0.01 0.0044 0]; %输入各环节参数 W0=zeros(9,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1; W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1; W(8,6)=1;W(9,7)=1; %输入系统连接矩阵 A=diag(P(:, 1));B=diag(P(:, 2));

C=diag(P(:, 3));D=diag(P(:, 4)); %生成系数矩阵 Q=B-D*W;Q1=inv(Q); R=C*W-A; V=C*W0;

A1=Q1*R; B1=Q1*V; %生成闭环系数矩阵 C1=[0 0 0 0 0 0 1 0 0];

u=1; %阶跃输入幅值

x0=zeros(9,1); t=0;y=0; %设置初值 x=x0;

for i=1:Tf/h

K1=A1*x+B1*u;

K2=A1*(x+h*K1/2)+B1*u; K3=A1*(x+h*K2/2)+B1*u; K4=A1*(x+h*K3)+B1*u;

x=x+h*(K1+2*K2+2*K3+K4)/6; y=[y;C1*x]; t=[t;t(i)+h]; end

plot(t,y)

仿真时间Tf=20 仿真时间Tf=5 计算步长h=0.015 计算步长h=0.015

404035353030252520201515101055002468101214161820000.511.522.533.544.55

仿真时间Tf=100 计算步长h=0.016 2x 103020-2-4-6-8-10-12-14-16-180204060801001206

（2）、 clc clear

Tf=input('仿真时间Tf='); %输入仿真时间 h=input('计算步长h='); %输入仿真步长 P=[1 0.01 0.1 0; 0.17 0.5 1 1.5; 1 0.01 1 0;

0 0.051 1 0.15; 1 0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0; 1 0.01 0.1 0;

1 0.01 0.0044 0]; %输入各环节参数 W0=zeros(9,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1; W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1;W(8,6)=1;W(9,7)=1; %输入系统连接矩阵 A=diag(P(:, 1));B=diag(P(:, 2));

C=diag(P(:, 3));D=diag(P(:, 4)); %生成系数矩阵 Q=B-D*W;Q1=inv(Q); R=C*W-A; V=C*W0;

A1=Q1*R; B1=Q1*V; %生成闭环系数矩阵 C1=[0 0 0 0 0 0 1 0 0];

u=1; %阶跃输入幅值

x0=zeros(9,1); t=0;y=0; %设置初值 x=x0;

for i=1:Tf/h

K1=A1*x+B1*u;

K2=A1*(x+h*K1/2)+B1*u; K3=A1*(x+h*K2/2)+B1*u; K4=A1*(x+h*K3)+B1*u;

x=x+h*(K1+2*K2+2*K3+K4)/6; y=[y;C1*x]; t=[t;t(i)+h]; end

plot(t,y)

7

仿真时间Tf=20 计算步长h=0.015

30252015105002468101214161820仿真时间Tf=5 计算步长h=0.015

30252015105000.511.522.533.544.55

SIMULINK仿真

8

0.2120.10.01s+10.17s+10.85sSaturation110.01s+10.15s+10.051s700.00167s+10.210.15s+1130s0.10.01s+1y1To Workspace0.00440.01s+1

9

0.2120.10.01s+10.17s+10.85ss+0.785s+0.1610.01s+110.15s+1s0.051s7010.2111300.00167s+10.5s+10.15s+10.1s+1s0.10.01s+1y1To Workspace0.00440.01s+1

实验二

1、已知某四阶非线性系统结构框图如图1所示。当系统输入幅值为10的阶跃信号时，试仿真分析系统输出响应：

1）在SIMULINK环境下建立仿真模型，输出仿真结果；

2）采用MATLAB语言编程实现仿真程序，与1)进行结果比较；

3）在第一环节前加饱和非线性环节，分析限幅C1=6、2时系统的输出响应，并与没限幅前比较其超调、峰值时间及调节时间。

解：（1）SIMULINK图形如下

10

s+0.785s+0.161s10.5s+110.1s+1s+0.785s+0.16Saturation1s10.5s+110.1s+1s+0.785s+0.16Saturation11s10.5s+110.1s+1

（2）

11

子程序：satura.m

function uc=satura(ur,s1) %ur=1;s1=5; if (ur>=s1) uc=s1;

elseif (ur<=-1*s1) uc=-1*s1; else uc=ur; end

子程序：dead.m

function uc=dead(ur,s1) %ur=-1.5; s1=0.5; if (abs(ur)>=s1) if (ur>0) uc=ur-s1; else uc=ur+s1; end

else uc=0; End 主程序： clc clear

Tf=20; %输入仿真时间 h=0.01; %输入仿真步长

z=input('非线性标志向量z='); %输入非线性标志向量

c1=input('饱和非线性环节限幅c1=');%输入饱和非线性环节参数 a1=6;%输入死区非线性环节参数 P=[0.16 1 0.785 1; 0 1 1 0; 1 0.5 1 0;

1 0.1 1 0]; %输入各环节参数

A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);

W0=zeros(4,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(4,4);W(1,4)=-1;W(2,1)=1;W(3,2)=1;W(4,3)=1; for i=1:4

if(A(i)==0); F(i)=1;

G(i)=h*C(i)/B(i);

H(i)=0.5*h*h*C(i)/B(i); if (D(i)==0)

C1(i)=1;D1(i)=0; else

C1(i)=1;D1(i)=D(i)/B(i); %求积分、比例积分环节离散系数

12

end else

F(i)=exp(-h*A(i)/B(i)); if (D(i)==0)

G(i)=(1- F(i))*C(i)/A(i);

H(i)=h*C(i)/A(i)-G(i)*B(i)/A(i); C1(i)=1;D1(i)=0; else

G(i)=(1- F(i))*D(i)/A(i);

H(i)=h*D(i)/A(i)-G(i)*B(i)/A(i); C1(i)=C(i)/D(i)-A(i)/B(i);

D1(i)=D(i)/B(i); %求惯性、比例惯性环节离散系数 end end end

?=diag(F(:, 1));G1=diag(G(:, 1));H1=diag(H(:, 1));?=diag(C1(:, 1));D2=diag(D1(:, 1));

%求各环节输入

Y=zeros(4,1);X=Y;y=zeros(1,4);r=10;Uk=zeros(4,1); t=0:h:Tf; N=length(t); %求采样点个数 for k=1:N-1

Ub=Uk; %保存前一次输入

Uk=W*Y+W0*r; %求当前个环节输入v for i=1:4;

if(z(i)~=0) if(z(i)==1)

Uk(i,1)=satura(Uk(i,1),c1); end

if(z(i)==2)

Uk(i,1)=dead(Uk(i,1),a1); end end end

Udot=(Uk-Ub)/h; %求当前各环节输入导数 Uf=2*Uk-Ub; %求下一拍输入 %求各环节状态变量及输出

X=F'.*X+G'.*Uk+H'.*Udot; Y=C1'.*X+D1'.*Uf; for i=1:4;

if(z(i)~=0) if(z(i)==3)

Y(i,1)=satura(Y(i,1),c1); end

if(z(i)==4)

13

Y(i,1)=dead(Y(i,1),a1); end end end

y=[y;Y']; end

y4=y(:,4); plot(t,y4)

非线性标志向量z=[0 0 0 0] 饱和非线性环节限幅c1=6

1614X: 2.74Y: 15.1712108642002468101214161820 （3）

非线性标志向量z=[1 0 0 0] 饱和非线性环节限幅c1=6

14

15X: 3.12Y: 14.7105002468101214161820

非线性标志向量z=[1 0 0 0] 饱和非线性环节限幅c1=2

14X: 4.91Y: 13.1812108642002468101214161820

2、 某控制系统如图所示，选择增益K的值，使系统阶跃响应的超调整量小于20%，且调节时间小于5 s。（教材P80习题2.11）

解：SIMULINK图形如下： K(s+1)(s+1)(s+8)(s+20)StepZero-Pole132s +3.2s +3.56sTransfer FcnScope 15

r(t)?2?1(t)1550(s+1)(s+1)5132s +3.2s +3.56sTransfer FcnSaturation2s+12s3 +2s +sTransfer FcnScope(s+8)(s+20)StepStepZero-PoleGainScope

670(s+1)(s+1)(s+8)(s+20)StepZero-Pole132s +3.2s +3.56sTransfer FcnScope

3、某非线性系统如图所示，试求 时系统的动态特性。 (P52

16

页例题2.7)

5 StepGain

解：SIMULINK图形如下：

MATLABFunctionMATLAB Fcn2s+12s3 +2s +sTransfer FcnScope1StepGainSaturation2s+132s +2s +sTransfer FcnScope

17

5ufcny2s+1s3 +2s2 +sStepGainTransfer FcnEmbeddedMATLAB Function52s+1s3 +2s2 +sStep1Gain1SaturationTransfer Fcn1其中：

function y = fcn(u) %#eml if u>=1 y=1;

elseif u<=-1 y=-1; else y=u; end

Scope

18

4、考虑如图所示的单位反馈系统，试用双线性变换和根匹配法对其进行仿真（r(t)=1(t),T=0.05)：

（1）通过“c2d”指令将系统各环节进行离散化；

（2）通过“tools—Control Design– Model Discretize”将连续系统离散化； （3）试绘制闭环频率特性曲线。

s+3.33101 57.662 s+24ss +10s+16Transfer Fcn3Transfer Fcn4Transfer Fcn1 解：（1） t=0.05;

n1=[1 3.33];d1=[1 24];[nd1,dd1]=c2dm(n1,d1,t,'tustin') n2=[10];d2=[1 0];[nd2,dd2]=c2dm(n2,d2,t,'tustin') n3=[1];d3=[1 10 16];[nd3,dd3]=c2dm(n3,d3,t,'tustin') 运行结果： nd1 =

0.6770 -0.5730

dd1 =

1.0000 -0.2500

nd2 =

0.2500 0.2500

dd2 =

1 -1

nd3 =

1.0e-003 *

0.4960 0.9921 0.4960

dd3 =

1.0000 -1.5714 0.6032

19

Step

（2）

57.66StepGaintustin(s+3.33)(s+24)Zero-Poletustin10sTransfer Fcn1tustins 2 +10s+16Transfer Fcn1Scope

20

（3）

57.66StepGain(s+3.33)(s+24)Zero-Pole10sTransfer Fcn12s +10s+16Transfer Fcn1Scope57.66Step1Gain10.6770-0.5730z -11.0000-0.2500z -1Discrete Filter0.2500+0.2500z -11-z -1Discrete Filter1-1-20.000496+0.0009921z +0.000496z -1-21.0000-1.5714z +0.6032z Discrete Filter2

21

实验三

1、（P89例3.1）某数字控制系统如图所示，其输入为单位阶跃信号，采样周期T=0.1s，试求系统在各采样时刻的输出y(t)。

r(t) e(t)Te(kT)11?0.5z?1数字控制器D(z)u(kT)T1?e?Tssu(t)1sx11s?1x2y(t)保持器Gh(s)被控对象G0(s)解：程序如下： clear all

m=2; n=3; % 明确输入、输出存储向量单元的维数； A=[2.4 -1.863 0.453]; % 差分方程输出项系数行向量 B=[0 0.005 0.005]; % 差分方程输入项系数行向量

R=zeros(m+1,1); % 建立参与递推运算的输入信号序列存储列向量

Y=zeros(n,1); % 建立参与递推运算的输出信号序列存储列向量 T=0.1; % 明确采样周期T=0.1s

Tf=15; M=Tf/T; % 设定仿真总时间为15s(共进行M=150次递推计算) yt=0;t=0; for k=1:M

R(k)=1; % kT采样时刻输入

R=[R(k);R(1:m)]; % 刷新参与递推运算的输入信号序列 yk=A*Y+B*R; % 递推运算

Y=[yk;Y(1:n-1)]; % 刷新参与递推运算的输出信号序列

yt=[yt,yk]; % yt为记载各采样(kT)时刻输出响应的行向量 t=[t,k*T]; % t为记载各采样(kT)时刻的行向量(与yt对应) end

plot(t,yt,'*m'); % 绘制各采样(kT)时刻的输出响应图 grid;

xlabel('time(s)'); ylabel('y(kT)');

1.51y(kT)0.5005time(s)1015

22

2、（P92例3.2）数字控制系统同例3.1，采样周期T=0.1s，采用双重循环方法求系统在单位阶跃信号r(t)=1(t) 作用下的状态响应设初始状态。设初始状态

解：程序如下：

F=[1 0;0.01 0.99];

x1(0)?x2(0)?0G=[0.01;0];c=[0 1]; % 连续被控对象的离散化状态方程的系数阵 TF=15; % 设定总仿真时间为15s

T=0.1;h=0.01; % 设定内循环步长h=0.01,而外循环则以采样周期T=0.1为步长 N1=round(TF/T); % 离散部分计算的次数(外循环)

N2=round(T/h); % 在一采样周期内连续被控对象计算次数(内循环) x=[0;0]; %为状态变量设置初始值 u=0;t=0; xt=x;

for k=1:N1 % 外循环 y=c*x; r(k)=1;

e(k)=r(k)-y;

u=0.5*u+e(k); % 外循环中递推求数字控制器在实际采样时刻的输出(步长为T) for j=1:N2 % 内循环

x=F*x+G*u; % 内循环中递推求受控对象状态响应(步长为h) xt=[xt,x]; % xt为记载各步状态响应的矩阵

t=[t,j*h+(k-1)*T]; % t为记载各采样(kT)时刻的行向量(与yt对应) end end

x1t=xt(1,:);x2t=xt(2,:); %输出状态变量x1,x2 plot(t,x1t,':k',t,x2t,'*m'); legend('x1', 'x2'); grid;

xlabel('time(s)');

21.81.61.41.210.80.60.40.20 05time(s)1015x1x2

23

3、采样控制系统数字仿真实验，系统采样周期为0.1s。求采用以下采样周期进行系统仿真，分析单位阶跃输入时的响应特性：1）T=0.1s， 2) T=0.04s

解：SIMULINK图形如下：

(z-0.98)(z-0.64)StepDiscreteZero-PoleZero-OrderHold1s+1Transfer Fcn1sIntegratorScope(z-0.98)(z-0.64)Step1DiscreteZero-Pole1Zero-OrderHold11s+1Transfer Fcn11sIntegrator1

24

1.% 这是例6.2的仿真程序 F=[1 0;0.01 0.99];

G=[0.01;0];c=[0 1]; % 连续被控对象的离散化状态方程的系数阵 TF=15; % 设定总仿真时间为15s

T=0.1;h=0.01; % 设定内循环步长h=0.01,而外循环则以采样周期T=0.1为步长

N1=round(TF/T); % 离散部分计算的次数(外循环)

N2=round(T/h); % 在一采样周期内连续被控对象计算次数(内循环) x=[0;0]; %为状态变量设置初始值 u=0;t=0; xt=x;

for k=1:N1 % 外循环 y=c*x; r(k)=1;

e(k)=r(k)-y;

u=0.5*u+e(k); % 外循环中递推求数字控制器在实际采样时刻的输出(步长为T)

for j=1:N2 % 内循环

x=F*x+G*u; % 内循环中递推求受控对象状态响应(步长为h) xt=[xt,x]; % xt为记载各步状态响应的矩阵

t=[t,j*h+(k-1)*T]; % t为记载各采样(kT)时刻的行向量(与yt对应) end end

x1t=xt(1,:);x2t=xt(2,:); %输出状态变量x1,x2 plot(t,x1t,':k',t,x2t,'*m'); legend('x1', 'x2'); grid;

xlabel('time(s)');

2.% 求连续被控对象离散化状态方程的程序(基于连续系统离散相似算法) h=0.01; % 虚拟采样周期为0.01s A=[0 0;1 -1]; B=[1;0];

[F,G]=c2d(A,B,h) % 输入端采用零阶保持器

3.% 这是例7.1的仿真程序 clear all

m=2; n=3; % 明确输入、输出存储向量单元的维数； A=[2.4 -1.863 0.453]; % 差分方程输出项系数行向量 B=[0 0.005 0.005]; % 差分方程输入项系数行向量

25

G(i)=(1- F(i))*D(i)/A(i);

H(i)=h*D(i)/A(i)-G(i)*B(i)/A(i); C1(i)=C(i)/D(i)-A(i)/B(i);

D1(i)=D(i)/B(i); %求惯性、比例惯性环节离散系数 end end end

?=diag(F(:, 1));G1=diag(G(:, 1));H1=diag(H(:, 1)); ?=diag(C1(:, 1));D2=diag(D1(:, 1));

%求各环节输入

Y=zeros(4,1);X=Y;y=zeros(1,4);r=10;Uk=zeros(4,1); t=0:h:Tf; N=length(t); t1=0; %求采样点个数 for k=1:N-1

Ub=Uk; %保存前一次输入

Uk=W*Y+W0*r; %求当前个环节输入v Udot=(Uk-Ub)/h; %求当前各环节输入导数 Uf=2*Uk-Ub; %求下一拍输入 %求各环节状态变量及输出

X=F'.*X+G'.*Uk+H'.*Udot; Y=C1'.*X+D1'.*Uf; y=[y;Y']; end

y4=y(:,4); plot(t,y4)

10. function uc=satura(ur,s1) %ur=1;s1=5; if (ur>=s1) uc=s1;

elseif (ur<=-1*s1) uc=-1*s1; else uc=ur; end

11.clc clear

Tf=input('仿真时间Tf='); %输入仿真时间 h=input('计算步长h='); %输入仿真步长

a=input('ASR参数a='); %输入转速调节器参数 P=[1 0.01 0.1 0; 0 0.85 1 a;

31

1 0.01 1 0;

0 0.051 1 0.15; 12.clc clear u=1;

x0=0; t=0;y=0; x0=[0;0;0];

A=[0 1 0;0 0 1;-22.06 -27 -10]; B=[0;0;1]; C=[40.6 0 0];

Tf=input('仿真时间Tf='); h=input('计算步长h='); x=x0;

for i=1:Tf/h K1=A*x+B*u;

K2=A*(x+h*K1/2)+B*u; K3=A*(x+h*K2/2)+B*u; K4=A*(x+h*K3)+B*u;

x=x+h*(K1+2*K2+2*K3+K4)/6; y=[y;C*x]; t=[t;t(i)+h]; end

plot(t,y)

1 0.0067 70 0; 1 0.15 0.21 0; 0 1 130 0; 1 0.01 0.1 0;

1 0.01 0.0044 0]; %输入各环节参数 W0=zeros(9,1);W0(1,1)=1; %输入外部链接矩阵

W=zeros(9,9);W(2,1)=1;W(2,9)=-1;W(3,2)=1;W(4,3)=1; W(4,8)=-1;W(5,4)=1;W(6,5)=1;W(6,7)=-0.212;W(7,6)=1;W(8,6)=1;W(9,7)=1; %输入系统连接矩阵 A=diag(P(:, 1));B=diag(P(:, 2));

C=diag(P(:, 3));D=diag(P(:, 4)); %生成系数矩阵 Q=B-D*W;Q1=inv(Q); R=C*W-A; V=C*W0;

A1=Q1*R; B1=Q1*V; %生成闭环系数矩阵 C1=[0 0 0 0 0 0 1 0 0];

u=1; %阶跃输入幅值

x0=zeros(9,1); t=0;y=0; %设置初值 x=x0;

for i=1:Tf/h

K1=A1*x+B1*u;

32

K2=A1*(x+h*K1/2)+B1*u; K3=A1*(x+h*K2/2)+B1*u; K4=A1*(x+h*K3)+B1*u;

x=x+h*(K1+2*K2+2*K3+K4)/6; y=[y;C1*x]; t=[t;t(i)+h]; end

plot(t,y)

13.clc clear u=1;

x0=[0 0 0]; C=[40.6 0 0]'; tspan=[0,10];

[t,x]=ode45('funexam2',tspan,x0); y=x*C;

plot(t,y);

14.clc clear

Tf=input('仿真时间Tf='); %输入仿真时间 h=input('计算步长h='); %输入仿真步长 x0=1/3; t=0; x=x0; %输入初始值 for i=1:Tf/h K1=-30*x0; x1=x0+h*K1;

x=[x;x1]; %变量以向量形式保存

t=[t;t(i)+h]; %对应时刻以向量形式保存 x0=x1; end

[t,x] %以数据形式输出 plot(t,x) %以曲线形式输出

15.clear

x0=[620,10,70]; % 置状态变量初值

tspan=[0,30]; % 置仿真时间区间

[t,x]=ode45('myf',tspan,x0); % 调用ode45求仿真解

plot(t,x(:,1),t,x(:,2),'--',t,x(:,3),':'); % 用不同的线型绘制仿真结果曲线

xlabel('time(天) t0=0, tf=30'); % 对t-x轴进行标注 ylabel('x(人):x1(0)=620,x2(0)=10;x3(0)=70');

33

legend('x1','x2','x3'); grid;

16.function dx=myf(t,x) dx = zeros(3,1);

dx(1)=-0.001*x(1)*x(2); % 第一个微分方程

dx(2)=0.001*x(1)*x(2)-0.072*x(2); % 第二个微分方程 dx(3)=0.072*x(2);

17.clear all

num=[0.005,0.005,0]; % 脉冲传递函数分子多项式按z的降幂系数 % 排列的行向量

den=[1,-2.4,1.863,-0.453]; % 脉冲传递函数分母多项式按z的降幂系数 % 排列的行向量

[yk,x,n]=dstep(num,den); % yk为存放输出离散序列的数组,n为dstep % 函数自动设定的采样点数 T=0.1; % 已知系统采样周期为0.1s for k=1:n

plot(k*T,yk(k),'*k'); % k为采样序列号,k*T为第k次采样对应的时刻 hold on end

xlabel('时间(s)'); grid hold off

18. function dy = vdp1000(t,y)

dy = zeros(2,1); % a column vector dy(1) = y(2);

dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);

34

本文来源：https://www.bwwdw.com/article/oajg.html