控制系统仿真技术（二）-连续系统的数字仿真二

太原理工大学控制系统仿真技术实验报告
连续系统的数字仿真
1.分别利用欧拉法和预估-校正法求下图所示系统的阶跃响应，并对其结果进行比较。

%欧拉法求阶跃响应
r=2;num0=8;den0=[1 3 0];
num1=1;den1=1;
[num,den]=feedback(num0,den0,num1,den1);
[A,B,C,D]=tf2ss(num,den);
Tf=input('仿真时间Tf=');h=input('计算步长h=');
x=[zeros(length(A),1)];y=0;t=0;
for i=1:Tf/h
    x=x+h*(A*x+B*r);
    y=[y;C*x];t=[t;t(i)+h];
end
plot(t,y);

%预估-校正法求阶跃响应
r=2;num0=8;den0=[1 3 0];
num1=1;den1=1;
[num,den]=feedback(num0,den0,num1,den1);
[A,B,C,D]=tf2ss(num,den);
Tf=input('仿真时间Tf=');h=input('计算步长h=');
x=[zeros(length(A),1)];y=0;t=0;
for i=1:Tf/h
    x1=x+h*(A*x+B*r);
    x=x+h*((A*x+B*r)+(A*x1+B*r))/2;
    y=[y;C*x];t=[t;t(i)+h];
end
plot(t,y);

%四阶龙格库塔法求阶跃响应
r=2;num0=8;den0=[1 3 0];
num1=1;den1=1;
[num,den]=feedback(num0,den0,num1,den1);
[A,B,C,D]=tf2ss(num,den);
Tf=input('仿真时间Tf=');h=input('计算步长h=');
x=[zeros(length(A),1)];y=0;t=0;
for i=1:Tf/h
    K1=A*x+B*r;
    K2=A*(x+h*K1/2)+B*r;
    K3=A*(x+h*K2/2)+B*r;
    K4=A*(x+h*K3)+B*r;
    x=x+h*(K1+2*K2+2*K3+K4)/6;
    y=[y;C*x];t=[t;t(i)+h];
end
plot(t,y);

2. 假设某一系统由四个典型环节组成，如下图所示：
分别用四阶龙格-库塔法和离散相似法求输出量 y 的动态响应曲线。
程序代码：（四阶—龙格库达法）

%针对系统结构图的四阶-龙格库塔法
r=10;P=[0.1  1 0.5 1;0 1 1 0;2 1 2 0;10 1 10 0];
W=[0 0 0 -1;1 0 0 0;0 1 0 0;0 0 1 0];
W0=[1;0;0;0];Wc=[0 0 0 1];
Tf=input('仿真时间Tf=');h=input('计算步长h=');
A1=diag(P(:,1));B1=diag(P(:,2));C1=diag(P(:,3));D1=diag(P(:,4));
H=B1-D1*W;Q=C1*W-A1;
A=inv(H)*Q;B=inv(H)*C1*W0;
x=[zeros(length(A),1)];y=[zeros(length(Wc(:,1)),1)];
 t=0;
 for i=1:Tf/h
     K1=A*x+B*r;
     K2=A*(x+h*K1/2)+B*r;
     K3=A*(x+h+K2/2)+B*r;
     K4=A*(x+h*K3)+B*r;
     x=x+h*(K1+2*K2+2*K3+K4)/6;
     y=[y;Wc*x];t=[t;t(i)+h];
 end
 plot(t,y);

（离散相似法）
即将下题目3中的P修改为
[0.1 1 0.5 1 0 0;0 1 1 0 0 0;2 1 2 0 0 0;10 1 10 0 0 0]

3.已知非线性系统如下图所示，试利用连续系统按环节离散化的数字仿真方法求 输出量 y 的动态响应，并与无非线性环节情况进行比较。

%ex2_23.m
%1前有饱和，2前有死区，3前有滞环，4是前有继电器，5后有饱和
%6后有死区，7后有滞环，8后有继电器
%注意P为从低到高
clear;
R=10;
P=[0.1 1 0.5 1 5 5;0 1 1 0 1 5;2 1 2 0 0 0;10 1 10 0 0 0];
W=[0 0 0 -1;1 0 0 0;0 1 0 0;0 0 1 0];
W0=[1;0;0;0];Wc=[0 0 0 1];
Tf=input('仿真时间Tf=');T=input('计算步长h=');
A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);FZ=P(:,5);S=P(:,6);
n=length(A);
for i=1:n
    if(A(i)~=0)
        if(B(i)==0)
            E(i)=0;F(i)=0;G(i)=0;H(i)=0;
            L(i)=(C(i)+D(i)/T)/A(i);Q(i)=-D(i)/(A(i)*T);
        else
            E(i)=exp(-A(i)*T/B(i));
            F(i)=(D(i)/B(i)-C(i)/A(i))*((1-E(i))*B(i)/(A(i)*T)-1);
            G(i)=(D(i)/B(i)-C(i)/A(i))*(1+(E(i)-1)*(1+B(i)/(A(i)*T)));
            H(i)=1;L(i)=D(i)/B(i);Q(i)=0;
        end
    else
        if(B(i)~=0)
            E(i)=1;F(i)=0.5*C(i)*T/B(i);G(i)=F(i);
            H(i)=1;L(i)=D(i)/B(i);Q(i)=0;
        else
            disp('A(i)=B(i)=0');
        end
    end
end
x=[zeros(length(A),1)];x0=x;z=x;
u=[zeros(length(A),1)];u0=u;
y=[zeros(length(Wc(:,1)),1)];t=0;
for j=1:Tf/T
    u1=u;u=W*x+W0*R;
    for i=1:n
        if(FZ(i)~=0)
            if(FZ(i)==1) 
                u(i)=saturation1(u(i),S(i));
            end
            if(FZ(i)==2) 
                u(i)=deadzone1(u(i),S(i));
            end
            if(FZ(i)==3)
                [u(i),u0(i)]=backlash1(u0(i),u(i),u1(i),S(i));
            end
            if(FZ(i)==4)
                u(i)=sign1(u(i),S(i));
            end
        end
    end
   x1=x;
   for i=1:n
       z(i)=E(i)*z(i)+F(i)*u(i)+G(i)*u1(i);
       x(i)=H(i)*z(i)+L(i)*u(i)+Q(i)*u1(i);
   end
    for i=1:n
        if(FZ(i)~=0)
            if(FZ(i)==5)
                x(i)=saturation1(x(i),S(i));
            end
            if(FZ(i)==6)
                x(i)=deadzone1(x(i),S(i));
            end
            if(FZ(i)==7)
                [x(i),x0(i)]=backlash1(x0(i),x(i),x1(i),S(i));
            end
            if(FZ(i)==8)
                x(i)=sign1(x(i),S(i));
            end
        end
    end
    y=[y;Wc*x];t=[t;t(j)+T];
end
plot(t,y)

注意这个程序里面还嵌套着saturation1，deadzone1，backlash1，sign1函数，所以还需要再相应创建这些函数

%饱和非线性
%saturation1.m
function x = saturation1(u,s)
if(abs(u)>=s)
    if(u>0) 
        x=s;
    else
        x=-s;
    end
else
    x=u;
end

%死区非线性
%deadzone1.m
function x = deadzone1(u,s)
if(abs(u)>=s)
    if(u>0)
        x=u-s;
    else
        x=u+s;
    end
else
    x=0;
end

%滞环非线性
%backlash1.m
function x = backlash1(u1,u,x1,s)
if(u>u1)
    if((u-s)>=x1)
        x=u-s;
    else
        x=x1;
    end
else if(u<u1)
        if((u+s)<=x1)
            x=u+s;
        else
            x=x1;
        end
    else
        x=x1;
    end
end
u1=u;

%继电器非线性
%sign1.m
function x = sign1(u,s)
if(u>0)
    x=s;
end
if(u<0)
    x=-s;
end