控制系统仿真技术（三）-采样系统的数字仿真

太原理工大学控制系统仿真技术实验报告
已知采样系统结构如下图所示。

试利用采样控制系统的数字仿真方法，当采样周期T=0.1s，且初始状态为零时，选取不同的计算步长求离散系统的单位阶跃响应。

%ex3_1.m
%采样系统的数字仿真
%其中Gr=[2.72 -1],Fr=[0.717]如果有离散过程2.72-z（-1）/1+0.717z(-1)
clear;
R=1;
Gr=[1];Fr=[0];
P=[1 1 1 0 0 0;2 1 1 0 0 0];
W=[0 0;1 0];
W0=[1;0];Wc=[0 1];
Tf=input('仿真时间Tf=');Tm=input('采样周期Tm=');T=input('计算步长h=');
A=P(:,1);B=P(:,2);C=P(:,3);D=P(:,4);FZ=P(:,5);S=P(:,6);
n=length(A);n1=length(Fr);m1=length(Gr);
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;
Ur=[zeros(n1,1)];Er=[zeros(m1,1)];
for ij=0:Tf/Tm
    e=R-x(n);Er=[e;Er(1:m1-1)];
    ur=-Fr*Ur+Gr*Er;Ur=[ur;Ur(1:n1-1)];
    for j=1:Tm/T
        u1=u;u=W*x+W0*ur;
         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(length(t))+T];
    end
end
plot(t,y)