自控MATLAB时域频域分析仿真

自动控制原理MATLAB仿真，求解传递函数、求解时域响应，串联校正法

串联校正法求法与教材（胡寿松）有所差异

%%%%%%%%%%%%%%%%%%%%%%%定义开环传递函数，求解闭环传递函数%%%%%%%%%%%%%%%%%%%%%%%%%
k = 5;
z = -20;
p = [0,-4.6,-1];
sys = zpk([z],[p],[k])      %以乘积形式呈现
[num,den] = zp2tf(z,p,k);       %转化位传递函数模型
g = tf(num,den);
sys = feedback(num1,den1,num2,den2,-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%求解时域响应%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num = [1.9691,5.0395];
den = [1,0.5572,0.6106];
step(num,den);      %单位脉冲响应
impulse(num,den);       %单位阶跃响应

%求解系统单位阶跃响应的性能指标程序
num=[22.68];
den=[1,9.2,7.36,5.376];
t=0:0.01:20;
step(num,den,t)         %单位阶跃响应
[y,x,t]=step(num,den,t)
maxy=max(y)         %峰值
yss=y(length(t))        %稳态值
pos=100*(maxy-yss)/yss      %超调量
for i=1:length(t)
    if y(i)==maxy
        n=i;end
end
tp=(n-1)*0.01       %峰值时间
y1=1.05*yss         %求解调节时间时的上限值
y2=0.95*yss         %求解调节时间时的下限值
i=2001              %从稳态往前递推
while i>0           %求解调节时间
    i = i-1;
    if y(i)>=y1|y(i)<=y2;m=i;break;
    end
end
ts = (m-1)/100;      %调节时间

%%%%%%%%%%%%%%%%%%%%%%频率分析法：串联超前校正%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%第一种情况：系统明确要求了校正后的截止频率
clear; close all; clc
G0=tf(10,[conv([1,0],[1,1])]);%待校正系统开环传函
bode(G0);%画出待校正系统的波特图
grid
hold on
margin(G0)%画出待校正系统的幅相特性曲线
[Gm_ori,Pm_ori,Wcg_ori,Wcp_ori]=margin(G0); %Gm为幅值裕度，Pm为相角裕度，Wcg为穿越频率，Wcp为截止频率

Wm = 4.8;%选取的截止频率
%根据20lg|G0(jWc)|=L=-10lg(a)求得参数a和L
a = 10^(-1*L/10);   %计算a值
T = 1/(Wm*sqrt(a)); %计算T值
Gc = tf([a*T,1],[T,1]);  %计算超前网络传递函数
G = G0*Gc;  %计算校正后的系统开环传递函数
margin(G)

%%%%%%%%%%%%%%%%%&第二种情况：系统只对稳态精度提出要求，没有对截止频率提出要求
clear; close all; clc
G0=tf(10,[conv([1,0],[1,1])]);%待校正系统开环传函
bode(G0);%画出待校正系统的波特图
grid
hold on
margin(G0)%画出待校正系统的幅相特性曲线
[Gm_ori,Pm_ori,Wcg_ori,Wcp_ori]=margin(G0); %Gm为幅值裕度，Pm为相角裕度，Wcg为穿越频率，Wcp为截止频率

theta = gama-gama_oral+10;%theta为最大超前相角、gama为题目要求的相角裕度、gama_oral为原系统相角裕度
a = (1+sin(theta))/(1-sin(theta));
%根据20lg|G0(jWc)|=-10lg(a)求得截止频率Wc
T = 1/(Wc*sqrt(a)); %计算T值
Gc = tf([a*T,1],[T,1]);  %计算超前网络传递函数
G = G0*Gc;  %计算校正后的系统开环传递函数
margin(G)

%%%%%%%%%%%%%%%%%%%%%%%频率分析法：串联滞后校正%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; close all; clc
G0=tf(10,[conv([1,0],[1,1])]);%待校正系统开环传函
bode(G0);%画出待校正系统的波特图
grid
hold on
margin(G0)%画出待校正系统的幅相特性曲线
[Gm_ori,Pm_ori,Wcg_ori,Wcp_ori]=margin(G0); %Gm为幅值裕度，Pm为相角裕度，Wcg为穿越频率，Wcp为截止频率

theta = -180+gama+10;%theta为最大超前相角、gama为题目要求的相角裕度
%根据Wc的相角等于theta，求出截止频率Wc
%根据20lg|G0(jWc)|+20lg(b)=0，求出参数b
T = 1/(0.1*Wc*b); %计算T值
Gc = tf([b*T,1],[T,1]);  %计算超前网络传递函数
G = G0*Gc;  %计算校正后的系统开环传递函数
margin(G)