《系统建模与仿真》实验代码（含题目）

LuanYu_

已于 2024-03-14 22:35:28 修改

阅读量720

点赞数 10

分类专栏：大学课程资料文章标签：课程设计笔记经验分享

于 2024-03-14 22:13:47 首次发布

本文链接：https://blog.csdn.net/Luan__Yu/article/details/136723280

版权

大学课程资料专栏收录该内容

13 篇文章 1 订阅

订阅专栏

实验五连续系统的仿真

（1）已知某系统的状态方程和输出方程为

$\dot{X(t)}=\begin{bmatrix} -8 & 1& 0\\ -19& 0& 1\\ -12& 0 & 0 \end{bmatrix}\cdot X(t)+\begin{bmatrix} 0\\ 4\\ 10 \end{bmatrix}\cdot u(t)$

$y(t)=\begin{bmatrix} 1 &0 &-3 \end{bmatrix}\cdot X(t)$

其中系统输入为 $u(t)=1(t)$ ，初始状态为 $X(0)=\begin{bmatrix} 0 &0 &0 \end{bmatrix}^{T}$ ，取步长h=0.01，试分别用欧拉法和四阶龙格-库塔法，求t=0.5 时y(0.5)的值。

（2）求上述线性系统的传递函数G(s)。对传递函数G(s)，应用双线性变换法，求出G(z)。假如采样周期T=0.01，且y(0)=0, y(0.01)=-0.3, y(0.02)=-0.6，求y(0.5)的值。

（1）

clear all
A=[-8 1 0;-19 0 1;-12 0 0];
B=[0;4;10];
C=[1 0 -3];
x=[0;0;0]; t=0;Y=0; x0=x; Y0=0;
h=0.01; T=0.5;
for i=1:T/h
    x0=x0+h*(A*x0+B);%应用欧拉法求解
    y0=C*x0;
    Y0=[Y0 y0];
    k1=A*x+B; %应用四阶龙格库塔法求解
    k2=A*(x+h/2*k1)+B;
    k3=A*(x+h/2*k2)+B;
    k4=A*(x+h*k3)+B;
    x=x+h/6*(k1+2*k2+2*k3+k4);
    y=C*x;
    Y=[Y y];
    t=[t i*h];
end
plot(t,Y0,'-',t,Y,'-')
y
y0

运行结果：

y =
  -13.1613
y0 =
  -13.2114

（2）

clear all
A=[-8 1 0; -19 0 1; -12 0 0];
B=[0; 4; 10];   C=[1 0 -3];
D=0;
x=[0; 0; 0];   t=0;Y=0; xo=x; Yo=0;
h=0.01;  T=0.5;
T=0.01;  %输入离散化步长
[num,den]=ss2tf(A,B,C,D,1);
sys=tf(num,den)  %转为传递函数的形式
sysd=c2d(sys,T,'tustin') %用双线性变换法将系统传递函数转换成脉冲传递函数
A=sysd.den{1};
B=sysd.num{1};
A=A(2:end);
L=length(B);
U=zeros(L,1);  Y=[-0.6;-0.3;0];
yt=0;t=0;
for k=3:0.5/T
    u(k)=1;
    U=[u(k); U(1:L-1)];  %刷新参与递推运算的输入信号序列
    y1=-A*Y+B*U;  %递推计算
    Y=[y1; Y(1:L-2)];  %刷新参与递推运算的输出信号序列
    yt=[yt; y1];
    t=[t; k*T];
end
plot(t,yt)
xlabel('t')
ylabel('y(t)')
grid on
y1

运行结果：

sys =
   -30 s^2 - 236 s - 416
  -----------------------
  s^3 + 8 s^2 + 19 s + 12
Continuous-time transfer function.
sysd =
 
  -0.1499 z^3 + 0.1383 z^2 + 0.1497 z - 0.1385
  --------------------------------------------
       z^3 - 2.921 z^2 + 2.844 z - 0.9231
采样时间: 0.01 seconds
离散时间传递函数。
y1 =
  -22.7220

实验六非线性系统的仿真

（1）描述Lorenz 系统的非线性微分方程为

$\left\{\begin{matrix} \dot{x}=10(y-x)\\ \dot{y}=28x-xz-y\\ \dot{z}=xy-\frac{8}{3}z \end{matrix}\right.$

已知初始状态为 $x(0) =y(0) =z(0) =0.1$ ，调用 ode45( )函数，试求该动力学模型在 $t\varepsilon [0, 20]$ 的数值解，并绘制变量 x(t), y(t), z(t)的三维空间曲线。

（2）搭建Simulink 仿真模型，绘制上述Lorenz 动力学模型中x(t), y(t), z(t) 的三维空间曲线。

（1）

function dx=sysexa(t,x)
  dx(1)=10*(x(2)-x(1));
  dx(2)=28*x(1)-x(1)*x(3)-x(2);
  dx(3)=x(1)*x(2)-(8/3)*x(3);
  dx=dx';
end

x0=[0.1 0.1 0.1];
tspan=[0,20];
[t,x]=ode45('sysexa',tspan,x0);
plot3(x(:,1),x(:,2),x(:,3))

运行结果：

（2）

plot3(out.simout,out.simout1,out.simout2)

运行结果：

实验七模型近似与 PID 控制器设计

（1）已知某被控对象的传递函数 $G_{0}(s)$ 为

$G_{0}(s)=\frac{12(s^2-3s+6)}{(s+1)(s+5)(s^2+3s+6)(s^2+s+2)}$

试对该高阶模型，求其带有时间延迟的一阶近似模型 $G_{1}(s)$ 。观察 $G_{0}(s)$ 和 $G_{1}(s)$ 的开环单位阶跃响应曲线，比较原模型 $G_{0}(s)$ 和近似模型 $G_{1}(s)$ 的接近程度。

（2）若将近似模型 $G_{1}(s)$ 作为开环传递函数，构成单位负反馈控制系统，绘出该闭环控制系统的单位阶跃响应，读出其超调量和调节时间；采用频域分析方法，求出系统的幅值裕度和相角裕度。

（3）利用教材中表4-1 的Z-N 整定公式，针对近似模型 $G_{1}(s)$ ，设计 PID 控制器 $D_{1}(s)$ 。

（4）使用Matlab 的pidtune()函数，针对近似模型 $G_{1}(s)$ ，设计 PID控制器 $D_{2}(s)$ 。分别从控制系统的超调量、调节时间、幅值裕度和相角裕度等方面，比较第3步和第4 步所设计控制系统的性能.

（5）使用Matlab 的交互式PID 控制器设计工具，为 $G_{1}(s)$ 设计 PID 控制器 $D_{3}(s)$ ，要求调节时间小于10s。记录此时的PID 控制器 $D_{3}(s)$ 的参数。

Lsbianshi函数

#Lsbianshi函数
function [gs,k,l,t]=lsbianshi(y,t)
fun=@(x,t)x(1)*(1-exp(-(t-x(2))/x(3))).*(t>x(2));
x=lsqcurvefit(fun,[1 2 3],t,y);
k=x(1); l=x(2); t=x(3)
gs=tf(k,[t 1],'iodelay',l);

（1）

s=tf('s');
sys=(12*(s^2-3*s+6))/((s+1)*(s+5)*(s^2+3*s+6)*(s^2+s+2));
[y,t]=step(sys);
gs=lsbianshi(y,t)
[y1,t1]=step(gs);
plot(t,y,'r')
hold on
plot(t1,y1,'k')

运行结果：

t =
    0.7241
gs =
                    1.217
  exp(-2.21*s) * ------------
                 0.7241 s + 1
Continuous-time transfer function.

（2）

figure(2)
gs1=feedback(gs,1);
[y2,t2]=step(gs1);
plot(t2,y2)
figure(3)
margin(gs);

（3）

k=1.217;l=2.21;t=0.7241;
a=k*l/t;
kp=1.2/a;ti=2*l;td=l/2;
ki=kp/ti;kd=kp*td;
numf=[kd kp ki]; denf=[1 0];
sysf=tf(numf,denf)
figure(4)
margin(sysf*gs)

运行结果：

sysf =
 
  0.357 s^2 + 0.3231 s + 0.07309
  ------------------------------
                s 
Continuous-time transfer function.

（4）

d=pidtune(gs,'pid');
figure(5)
step(feedback(gs*d,1),'k')
hold on
step(feedback(gs*sysf,1),'r')
figure(6)
margin(d*gs)

运行结果：

d =
             1 
  Kp + Ki * ---
             s 
  且 Kp = 0.348, Ki = 0.25
并联型的连续时间 PI 控制器。

（5）

m=pidtool(gs,'pid')

运行结果：

《系统建模与仿真》实验代码（含题目）

实验五 连续系统的仿真

实验六 非线性系统的仿真

实验七 模型近似与 PID 控制器设计

实验五连续系统的仿真

实验六非线性系统的仿真

实验七模型近似与 PID 控制器设计