自动控制原理线性系统的校正系统中, 计算串联校正装置的时候非常复杂..手算简直想死.

所以用matlab做了点微小的贡献…谢谢大家,code is here

使用方法

封装成了单个函数, 只要输入原始传函, 目标相角/幅角裕度, 返回校正后(可能多级校正)的传函, 以及bode图的前后对比

sys = tf([40],conv(conv([1,0],[0.2,1]),[0.0625,1]));
adjust(sys,30,15);

原始传函$$\frac{40}{0.0125 s^3 + 0.2625 s^2 + s}$$

校正装置$$\frac{0.009699 s^2 + 0.2209 s + 1}{0.0004482 s^2 + 0.0439 s + 1}$$

校正后的传函
$$\frac{0.388 s^2 + 8.836 s + 40}{5.603\cdot10^{-06} s^5 + 0.0006664 s^4 + 0.02447 s^3 + 0.3064 s^2 + s}$$

可以看到校正后裕度满足30,15的裕度条件

程序结构

四个函数(写的比较丑, 只是随写随用的, 所以冗余可能比较多

• 主函数 adjust

function csys = adjust(sys,w,h)

% test sys
% sys = tf([40],conv(conv([1,0],[0.2,1]),[0.0625,1]))
origin = surfeit(sys);%计算裕度
phase_forward = w - origin(1);%超前角

if phase_forward < 60 % 当超前角小于60时使用超前较正
    csys = forward(sys,w,h);
else % 否则使用滞后矫正
    csys = backward(sys,w,h);
end

[n,d] = tfdata(sys);
n = n{1};
n = n(find(n~=0));
d = d{1};
d = d(find(d~=0));

[n1,d1] = tfdata(csys(3));
n1 = n1{1};
n1 = n1(find(n1~=0));
d1 = d1{1};
d1 = d1(find(d1~=0));

csys(4) = tf(deconv(n1,n),deconv(d1,d));

bode(sys,csys(3))
legend('sys-origin','sys-result')
end

• 超前校正 forward

function csys = forward(sys,w,h)
w1 = 0.1:0.1:1000;
[mag,phase] = bode(sys,w1);
mag = 20*log10(mag);
mag = mag(:)';
f_mag = @(x) interp1(w1,mag,x);  %构造幅相函数
phase_plus = 10;
chara = surfeit(sys);
phase_chase = w - chara(1) + phase_plus;
sin_chase = sin(phase_chase * pi /180);
a = (1+sin_chase)/(1-sin_chase);
wc = fzero(@(x) f_mag(x)+10*log10(a), 15);
T = 1/(sqrt(a)*wc);
resc = tf([a*T,1],[T,1]) ;
res = resc * sys;
chara2 = surfeit(res);%矫正装置传递函数计算
if chara2(1) < w || chara2(2) < h
    an = forward(res,w,h);%不满足要求继续迭代, 多级串联
    csys(1) = an(1);
    csys(2) = an(2);
    csys(3) = an(3);
else
    csys(1) = sys;%满足直接返回结果
    csys(2) = resc;
    csys(3) = res;
end
end

• 滞后校正 bakcward

function csys = backward(sys,w,h)

res = zeros(1,2);
w1 = 0.1:0.1:1000;
[mag,phase] = bode(sys,w1);
mag = 20*log10(mag);
mag = mag(:)';
phase = phase(:)';
phase = 180 + phase;
f_mag = @(x) interp1(w1,mag,x);
f_phase = @(x) interp1(w1,phase,x); %构造幅相函数


chara = surfeit(sys);
wc = fzero(@(x) f_phase(x)-(w+9), 10);
b =fzero(@(x) f_mag(wc) + 20*log10(x), 0.1);
T = 1/(0.1*b*wc);

resc = tf([b*T,1],[T,1]) ;
res = resc * sys;
chara2 = surfeit(res); %矫正装置传递函数计算

if chara2(1) < w || chara2(2) < h
    an = backward(res,w,h);%不满足要求继续迭代, 多级串联
    csys(1) = an(1);
    csys(2) = an(2);
    csys(3) = an(3);
else
    csys(1) = sys;%满足直接返回结果
    csys(2) = resc;
    csys(3) = res;
end
end

• 裕度函数 surfeit

function res = surfeit(sys)
%计算传递函数幅值裕度, 相角裕度, 截止频率
res = zeros(1,2);
w = 0.1:0.1:1000;
[mag,phase] = bode(sys,w);
mag = 20*log10(mag);
mag = mag(:)';
phase = phase(:)';
phase = 180 + phase;
f_mag = @(x) interp1(w,mag,x);
f_phase = @(x) interp1(w,phase,x);
wc = fzero(f_mag, 15);
wx = fzero(f_phase , 15);
res(1) = f_phase(wc);
res(2) = -f_mag(wx);
res(3) = wc;
res(4) = wx;

end