autodesign.m

频率法校正程序

function index=autodesign(num,den,wm,kc,op)
%==========================================================================
%自动控制原理课程设计实验程序。
%输入变量num,den分别为开环环传递函数展开式中分子与分母系数
%校正采用频率法超前校正，wm为期望截止频率，当不实行校正时wm置0
%kc为根轨迹增益，kc取不同的值可得到不同增益校正后的传递函数，若不用则置1
%op为操作项，op为0时，画闭环阶跃响应曲线，op为1时，画波德图，op为2时画根轨迹，op为3时对原系统进行校正，且画校正后阶跃响应曲线，此时wm不能为0
%输出为时域和频域能指标，格式为index=[td tr tp ts siagama r h wc wx]
%当输出index矩阵有两行时，第一行为校正前的指标，等二行为校正后的指标
%当op不为3时，则不进行校正操作，index只有一行，为原系统性能指标
%==========================================================================
%==========================================================================
if op~=3
    wm=0;
end
go=tf(num,den);%求开环传递函数
go=kc*go;
[index,mag,w]=ide(go); %求当前系统各项时域和频域指标
%画阶跃响应曲线


%频率法校正
if wm~=0
    [wmmin,loc]=min(abs(w-wm));
    wm=w(loc);
    a=(1/mag(loc))^2;
    T=1/(wm*sqrt(a));
    gc=tf([a*T 1],[T 1]);
    gcgo=gc*go;
    index=[index;ide(gcgo)];%求校正后系统的各项指标
end

switch op
    case 0        
        go=go/(1+go);
        [y,t]=step(go,[0:0.1:40]);  %从0s,到100s,以0.1s为精度做阶跃响应实验
        plot(y);hold on;plot([0:400],ones(1,401),':r'); %画阶跃响应曲线
    case 1
        bode(go);grid on;%画波得图
    case 2
        rlocus(go); %做根轨迹
    case 3
        gcgo=gcgo/(1+gcgo);
        [y,t]=step(gcgo,[0:0.1:40]);  %从0s,到100s,以0.1s为精度做阶跃响应实验
        plot(y);hold on;plot([0:400],ones(1,401),':r'); %画阶跃响应曲线
end
        
function [index,mag,w]=ide(sys)%求该传递函数性能指标（频域及时域）
sysp=sys/(1+sys);
[y,t]=step(sysp,[0:0.1:100]);  %从0s,到100s,以0.1s为精度做阶跃响应实验             
len=length(y);
yend=y(len);                %取终值
ymax=max(y);                %取最大值（峰值）
sigama=(ymax-yend)/yend;    %计算超调量

%计算调节时间---------------------------------------
for ii=len:-1:1
    if abs(y(ii)-yend)>0.05*yend   
        ts=t(ii);               
        break
    end
end

for ii=2:len                      %计算峰值时间
    if y(ii)-y(ii-1)<=0
        tp=t(ii);
        break
    end
end

%计算上升时间-----------------------------------------
for ii=2:len           
    if y(ii)>=0.1*yend & y(ii-1)<0.1*yend
        tr1=t(ii);
    end
    if y(ii)>0.9*yend & y(ii-1)<=0.9*yend
        tr2=t(ii-1);
        break
    end
end
tr=tr2-tr1;                              
%计算延迟时间----------------------------------------
for ii=2:len
    if y(ii)>=0.5*yend & y(ii-1)<0.5*yend  
        td=t(ii); 
        break
    end
end
%频域分析--------------------------------------------         
[mag phase w]=bode(sys);
[magmin,loc]=min(abs(mag-1));
wc=w(loc);            %计算截止频率
r=180+phase(loc);     %计算相角裕度
[phasemin,loc]=min(abs(phase+180));
wx=w(loc);           %计算穿越频率
h=1/abs(mag(loc));  %计算幅值裕度
index=[td tr tp ts sigama r h wc wx];