Matlab做二阶函数

关于解答：

>> % Expert PID Control
% 二阶传递函数的阶跃响应
% 位置式
clc;
clear all;
close all;

n=500;                              % 设置离散点的个数 
Ts=0.001;                           % 设置离散的采样时间
epsilon=0.001;                      % 设置加入积分换进的精度ε
% 对连续时间的传递函数进行离散化
% 两种方式可以同样设置传递函数（连续时间）
% s=tf('s');
% Gs=(100+学号后两位)/(s^2+(10+学号后两位)*s);
Gs=tf(100+学号后两位,[1,(10+学号后两位),0]);
Gz=c2d(Gs,Ts,'z');
[num,den]=tfdata(Gz,'v');
time=zeros(1,n );                   % 离散时间
r=ones(1,n );                       % 离散输入R(s) r(t)
u1=zeros(1,n );                     % 离散专家PID控制器输出U(s) u(t)
u2=zeros(1,n );                     % 离散普通PID控制器输出U(s) u(t)
y1=zeros(1,n );                     % 专家PID控制下的系统输出
y2=zeros(1,n );                     % 普通PID控制下的系统输出
error1=zeros(1,n);                  % 专家PID控制下的系统误差
error2=zeros(1,n);                  % 普通PID控制下的系统误差
% 根据z变换后的公式可知需要
u1_1=0;u1_2=0;                      % u(k-1) u(k-2)
y1_1=0;y1_2=0;                      % y(k-1) y(k-2)
u2_1=0;u2_2=0;                      % u(k-1) u(k-2)
y2_1=0;y2_2=0;                      % y(k-1) y(k-2)
x1=[0,0,0]';                        % 位置式PID：x(1)=error(k);x(2)=(error(k)-error(k-1))/Ts;x(3)=Ts*Σen (n=1~k)
x2=[0,0,0]'; 
x1_2_1=0;
x2_2_1=0;
kp1=2.0;ki1=0.20;kd1=1.00;          % 专家PID原有PID参数
% kp1=3.5;ki1=0.05;kd1=0.25;          % 专家PID原有PID参数
kp2=2.0;ki2=0.20;kd2=1.00;          % 专家PID原有PID参数
% kp2=3.5;ki2=0.05;kd2=0.25;          % 普通PID原有PID参数
% kp2=3.0;ki2=0.05;kd2=0.1;         % 普通PID原有PID参数
error1_1=0;                         % 误差的第k个极值
error2_1=0;                         % 误差的第k个极值

% k1=5.0;                             % 增益放大系数
k1=2.0;  
k2=0.6;                             % 抑制系数

% 循环更新数值(专家PID)
for k=1:1:n  
    time(k)=k*Ts;                   %离散的时间
    u1(k)=kp1*x1(1)+kd1*x1(2)+ki1*x1(3);   %PID控制器
    
    %专家控制的规则     
    %规则1
    %当误差的
    if abs(x1(1))>0.80
        u1(k)=(3/4)*kp1;
    elseif abs(x1(1))>0.40
        u1(k)=(2/3)*kp1;
    elseif abs(x1(1))>0.20
        u1(k)=(1/5)*kp1;
    elseif abs(x1(1))>0.01
        u1(k)=(1/10)*kp1;
    end
    
    %规则2
    if (x1(1)*x1(2)>0)||(x1(2)==0)
        if abs(x1(1))>=0.05
           u1(k)=k1*(kp1*x1(1)+kd1*x1(2)+ki1*x1(3)); 
        else
           u1(k)=k2*(kp1*x1(1)+kd1*x1(2)+ki1*x1(3));
        end
    end
    
    %规则3
    if ((x1(1)*x1(2)<0)&&(x1(2)*x1_2_1>0))||(x1(1)==0)
        u1(k)=u1(k);
    end
    
    %规则4
    if (x1(1)*x1(2)<0)&&(x1(2)*x1_2_1<0)
        if abs(x1(1))>=0.05
            u1(k)=k1*kp1*error1_1;
        else
            u1(k)=k2*kp1*error1_1;
        end
    end
    
    %规则5
    if abs(x1(1))<=epsilon
        u1(k)=((5/6)*kp1)*x1(1)+((1/3)*ki1)*x1(3);
%         u1(k)=0.5*x1(1)+0.01*x1(3);
    end
    
    % 限制控制器的范围
    if u1(k)>=10
        u1(k)=10;
    end
    if u1(k)<=-10
        u1(k)=-10;
    end
    
    %线性化的传递函数
    y1(k)=-den(2)*y1_1-den(3)*y1_2+num(1)*u1(k)+num(2)*u1_1+num(3)*u1_2;
    error1(k)=r(k)-y1(k);
    
    %更新参数
    u1_2=u1_1;u1_1=u1(k);
    y1_2=y1_1;y1_1=y1(k);
    x1(1)=error1(k);
    x1_2_1=x1(2);
    x1(2)=(error1(k)-error1_1)/Ts;
    x1(3)=x1(3)+error1(k)*Ts;
    error1_1=error1(k);
end

% 循环更新数值(普通PID)
for k=1:1:n  
    time(k)=k*Ts;                   %离散的时间
    u2(k)=kp2*x2(1)+kd2*x2(2)+ki2*x2(3);   %PID控制器
    
    %线性化的传递函数
    y2(k)=-den(2)*y2_1-den(3)*y2_2+num(1)*u2(k)+num(2)*u2_1+num(3)*u2_2;
    error2(k)=r(k)-y2(k);
    
    %更新参数
    u2_2=u2_1;u2_1=u2(k);
    y2_2=y2_1;y2_1=y2(k);
    x2(1)=error2(k);
    x2_2_1=x2(2);
    x2(2)=(error2(k)-error2_1)/Ts;
    x2(3)=x2(3)+error2(k)*Ts;
    error2_1=error2(k);
end

figure(1);
plot(time,r,'b',time,y1,'r');
xlabel('time(s)');
ylabel('r,y');
legend('输入','专家PID输出');
figure(2);
plot(time,r-y1,'r');
xlabel('time(s)');
ylabel('error');
legend('专家PID误差');
figure(3);
plot(time,r,'b',time,y1,'r',time,y2,'k');
xlabel('time(s)');
ylabel('r,y1,y2');
legend('输入','专家PID输出','普通PID输出');