一些控制方法在matlab实现的步骤及代码

目录

一.物理建模得被控对象

1.传递函数

2.状态空间方程

做法一:直接构建ABCD矩阵

做法二 :由传递函数得到状态方程

二.系统分析

1.设计要求

2.在阶跃响应下开环传递函数输出波形

做法一:带figure

做法二:直接输出波形图

做法三:线性系统分析器

3.看看开环传递函数零极点

做法一:pole+zero指令

做法二:pzmap指令

4.得到闭环传递函数

5.对输出Y(s)进行因式分解以求得拉普拉斯反变换   

(若无特殊要求,此步省略)

三.控制方法与分析

1.根轨迹法

先找补偿器增益Kp,看看单纯用Kp能不能满足要求

超前补偿器

滞后补偿器

2.频率分析法

先看开环伯德图

设计比例控制器Kp

超前补偿器

滞后补偿器

3.状态空间法

观察能控性(全状态反馈控制器)

观察能观性(状态观测器)

设计全状态反馈控制器

预补偿器Nbar

4.PID

5.数字控制器

(1)做法一:状态方程+LQR 

(2)做法二:传递函数+根轨迹

(3)做法三:传递函数+PID+根轨迹

四.附加资料


主要包含根轨迹+频率分析+状态空间+数字控制器+PID。

拿到物理模型→系统分析→列出物理运动方程式→拉普拉斯变换→传递函数+状态空间方程→引入控制器→系统输出响应波形即可达到期望要求

其方法及步骤内容主要来自:湖南工业大学自动控制设计实验室

一.物理建模得被控对象

1.传递函数

%传递函数
s=tf('s');
P_pitch=(1.151*s+0.1774)/(s^3+0.739*s^2+0.921*s);

2.状态空间方程

做法一:直接构建ABCD矩阵

%状态空间方程
A = [-0.313 56.7 0;-0.0139 -0.426 0;0 56.7 0];
B = [0.232; 0.0203;0];
C = [0 0 1];
D = [0];
pitch_ss=ss(A,B,C,D);

做法二 :由传递函数得到状态方程

(1)自己写出来传递函数分子 分母

%自己写传递函数分子分母
num=[1.151,0.1774];
den=[1,0.739,0.921,0];
[A,B,C,D]=tf2ss(num,den) %num是开环传递函数分子 den是分母
pitch_ss_1=ss(A,B,C,D)

(2)直接由传递函数变来

pitch_ss_2=ss(P_pitch)%其中P_pitch是开环传递函数

注意:此命令得到状态方程的ABCD并不会直接保存在工作区 


二.系统分析

1.设计要求

  •  超调量
  • 上升时间
  • 稳定时间
  • 稳态误差

2.在阶跃响应下开环传递函数输出波形

做法一:带figure

%看看开环传递函数输出响应波形
t=[0:0.01:10];     %时间从0→10s,中间每间隔0.01就记录一下数值大小,即采样
figure(1);         %告诉软件,我要画一个图,先把图方框准备好
step(P_pitch,t);   %阶跃响应step下,开环传递函数的输出响应波形。
grid on;           %显示网格
xlabel('时间');   %横  纵坐标的名字
ylabel('飞机的俯仰角');
title('开环阶跃响应曲线');  %标题

做法二:直接输出波形图

step(P_pitch),grid  %不加;号直接显示输出
axis([0 10 0 5])%限制坐标范围。横坐标范围0-10  纵坐标范围0-5
xlabel('时间')
ylabel('飞机的俯仰角');
title('开环阶跃响应曲线');

做法三:线性系统分析器

%启动线性系统分析器,输入波形是step,研究对象是开环传递函数,时间为0→10.
linearSystemAnalyzer('step',P_pitch,0:0.01:10);
%可在分析器里看到输出阶跃响应波形

注意:线性系统分析器功能强大,例如

  • 可以更改输入响应。右键--线性仿真--设计信号。(或者在命令语句将step改为正弦波等等)
  • 右键可以更改想看的波形,可以改看零极点分布图,bode图,rlocus根轨迹图等等
  • %可以同时输入多个传递函数,观察波形分析
    rP_pitch = 0.1/(0.5*s+1) %随便设的
    linearSystemAnalyzer('step', rP_pitch, P_pitch,0:0.1:10);

3.看看开环传递函数零极点

做法一: pole zero指令

fprintf('开环传递函数极点为');
pole(P_pitch)    %开环极点 注意没有;号
zero(P_pitch)    %开环零点

做法二:pzmap指令

pzmap(P_pitch)   %同时看零极点
axis([-1 1 -1 1])

4.得到闭环传递函数

fprintf('闭环传递函数为');
sys_cl=feedback(P_pitch,1)


 5.对输出Y(s)进行因式分解以求得拉普拉斯反变换   

(若无特殊要求,此步省略)

输出Y(s)=传递函数G(s)×输入R(s)。此处输入R(s)=  1/s。

R=1/s;
Y=zpk(sys_cl*R)

                                        1.151 (s+0.1541)
①即Y(s)=            ------------------------------------------
                       s  (s+0.08805) (s^2 + 0.6509s + 2.015)

②当前目标是因式分解成如下形式

\frac{A}{s}+\frac{B}{s+0.08805}+\frac{Cs+D}{s^2 + 0.6509s + 2.015}

③现在求ABCD(也可以手动计算),若使用matlab计算,相关命令为residue

[r,p,k]=residue([1.151 0.1774],[1 0.739 2.072 0.1774 0])

得到             -0.2802 + 0.0800i                                -0.3255 + 1.3816i
                    -0.2802 - 0.0800i                                 -0.3255 - 1.3816i 
        r=         -0.4395 + 0.0000i                        p=    -0.0881 + 0.0000i                        k=[]
                    1.0000 + 0.0000i                                   0.0000 + 0.0000i
④因此A=1,B=-0.0881,现在求C D

[num,den]=residue(r(1:2),p(1:2),k);
tf(num,den)

⑤可得                                                                           

                -0.5605 s - 0.4036
                --------------------------
              s^2 + 0.6509 s + 2.015      

⑥因此C=-0.5605     D=  -0.4036            A B C D 全部得到,可以进行拉普拉斯反变换求解

syms s

A=1;                  
B=-0.0881;
C=-0.5605;
D=-0.4036;
F=A/s+B/(s+0.08805)+(C*s+D)/(s^2 + 0.6509*s + 2.015);

ilaplace(F)

⑦这样y(t)就得到了



三.控制方法与分析

1.根轨迹法

先找补偿器增益Kp,看看单纯用Kp能不能满足要求

(1)做法一:

%打开控制系统设计器
controlSystemDesigner('rlocus',P_pitch)
%命令中,VIEWS 必须为以下一个或多个项: 'rlocus'、'bode'、'nichols' 和 'filter'

 注意:控制系统设计器功能强大,例如

  • 重要:可同时查看根轨迹与输出响应,并且根轨迹的零极点可拖动,移动根轨迹上的零极点从而改变K,使得输出响应实时变化。
  • 针对不同的回路结构,可以在编辑架构里选择不同的回路
  • 点击控制器G,可实时查看当前回路增益Kp与开环传递函数
  • 重要:根轨迹中,右键--设计需求--新建,可添加要求。其中上升时间与固有频率有响应转换关系,要求添加完毕后,非阴影区(白色区)即是期望区域。
  • 输出响应也可以右键添加指标,显示当前的稳定值 上升时间等。

(2)做法二:

rlocus(P_pitch)
axis([-0.6 0 -0.6 0.6]);%坐标范围
sgrid(0.6,0.9)%第一个是阻尼比,第二个是固有频率(有公式算出来的)

根轨迹图有了,接下来找Kp。找Kp就必须先在根轨迹上选一个点。

[Kp,pole]=rlocfind(P_pitch)  %选择极点后对应的K值

命令执行后会让你在根轨迹上射一个点,选择此点后会显示①此点对应的Kp值+②离这个点最近的、实轴上的 极点坐标

以上两个方法都得到了Kp值,把这个Kp值乘进开环传递函数看看输出响应波形,如果单纯的Kp无能为力,则还有超前/滞后补偿器可以抢救。

超前补偿器

适用条件:当前的根轨迹与所要求的非阴影区(期望区域)(要求区域)不沾边,具体来说当前根轨迹范围完全在期望区域的右侧,就不得不向左移动根轨迹。(换句话说,当我们希望闭环的根往左移动(更稳定时),引入超前补偿即可)

传递函数

C(s)=\frac{s+z}{s+p}

(1)做法一:

继续在控制系统设计器中调节

  • 点击右上角预设项--选项--零点/极点/增益--确定
  • 右击空白处--编辑器补偿器--右击--添加零极点--超前--修改实极点 实零点的值(选择有技巧)(此例的零点选择为-0.9  极点选择为-3)

随后新根轨迹图就有了,拖动极点以达到期望的输出响应波形。

(2)做法二:

%带超前补偿器的根轨迹图
z0=0.9;
p0=3;
s=tf('s');
C_lag=(s+z0)/(s+p0);  %构建超前补偿器
rlocus(C_lag*P_cruise) %看根轨迹图
axis([-0.6 0 -0.6 0.6]);
sgrid(0.6,0.36)
%重新选Kp
[Kp,pole]=rlocfind(C_lag*P_cruise)
再用此Kp去求输出波形
Kp =200;
sys_cl=feedback(Kp*C_lag*P_cruise,1);
t=0:0.1:20;
step(r*sys_cl,t)

(3)做法三:

可以尝试不同类型的动态补偿器(添加零点和极点)

滞后补偿器

适用条件:响应速度都挺好,就是稳态误差太大,稳定值跟期望值差老远。因此滞后补偿器可以在减小稳态误差的同时,仍能满足暂态要求(不会过大的改变系统原有的特性)。

传递函数

C(s)=\frac{s+z}{s+p}

完全同超前补偿器,唯独不同之处在于z值与p值的选取,超前补偿器|z|<|p|,滞后补偿器|z|>|p|

例如 超前补偿器:                                                        滞后补偿器:

                              C(s)=\frac{s+3}{s+0.3}                                                   C(s)=\frac{s+0.3}{s+3}

其他步骤一模一样。


特殊情况:

在某些特殊情况下,例如什么补偿器都加了,共轭极点已经拉到 期望范围 的边缘了,但是输出波形还是不满足要求,有可能是因为实轴上的、靠近虚轴的第三极点作为主导极点,不禁减缓了系统响应,而且压制了共轭极点的影响,这个情况下,可以将共轭极点拉到 期望范围 外,依然会满足要求。


2.频率分析法

先看开环伯德图

(1)做法一:

%在控制系统设计器里看开环传函伯德图
controlSystemDesigner('bode', P_pitch)

(2)做法二: 

bode(P_pitch),grid

(3)做法三: 

margin(P_pitch),grid       %还可以看到增益裕度与相位裕度
%若增益裕度与相位裕度都为正,则系统稳定

设计比例控制器Kp

如果当前bode图与输出波形不符合要求,则有以下两种方法找到增益Kp

(1)做法一:比较适用于一阶系统

①稳态误差公式          error=\frac{1}{1+M_{w\rightarrow 0}}*100%

②在bode图上找到w→0时对应的幅值dB(例如-34dB),带入公式计算得到稳态误差,同时也可以再看看闭环输出响应来验证一下(Kp取1即可)

sys_cl=feedback(Kp*P_pitch,1);
step(sys_cl)

③减小稳态误差,即抬高稳态值,例如要求稳态误差小于2%,带入公式 可以计算得出期望幅值=33.8dB,因此从原本的-34dB→期望的33.8dB,一次性抬高67.8dB。再根据幅值与比例增益Kp之间的关系,可得Kp=2455.

④随后便可以将Kp=2455乘入传递函数,观察开传函bode图以及闭环输出波形

⑤特殊情况,由于这种做法Kp往往取的很大,在某些特殊控制情况下(例如汽车巡航控制,输出速度响应波形变化太快,现实里根本无法实现)(例如Kp过大导致系统有振荡,相位裕度太小,系统不稳定等),就不得不降低Kp,可以引入滞后补偿器。

(2)做法二:比较适用于多阶系统

①直接考虑相位裕度,一般为了得到满意的性能,相位裕度一般为30°~60°,这里直接拉满令Pm=60deg。180°-相位120°=相位裕度60°,因此直接找到相位120°所对应的点

②在原始开环bode图中找相位120°所对应的频率(例如=9.17rad/s),将此频率作为新的截止频率

③查看频率为9.17所对应的mag(幅值图的横坐标)与相角phase

[mag,phase,w] = bode(P_pitch,9.17)

④得到数据(举例)mag=0.0157   再将mag换算成幅值

20*log10(0.0157)

⑤得到幅值=-36.0820dB,为了使当前选中频率9.17rad/s为截止频率,就必须将此幅值-36.0820dB→0dB,因此需要抬高36.0820dB,换算成Kp值≈63.

⑥利用此Kp值再去看新bode图(Pm≈60)与输出阶跃响应。如果输出波形不满意,引入补偿器

超前补偿器

特别注意:使用补偿器的时候要从系统的开环传递函数从头设计,以上得到的Kp值全部扔掉。

使用条件:加快截止频率且增加相位裕度。它为系统增加了正相位,这样可以增加相位裕度,从而增加阻尼。与此同时,它还可以在高频率下增加开环频率响应的幅度,进而增强截止频率和整体速度(稳定时间降低)。

传递函数:

                        C(s)=K\frac{Ts+1}{\alpha Ts+1}        (\alpha <1)

步骤:

①先只引入一个K,先初步决定K=10,看看根轨迹+输出波形

K=10;
margin(K*P_pitch),grid
figure;
sys_cl=feedback(K*P_pitch,1);
step(sys_cl),grid

②在本例中,增加K=10会使得输出有振荡,这是因为相位裕度变小了(本例Pm=10.4deg),不过没关系,超前补偿器可以增大阻尼、消除振荡。

③现在选择α,α需要由最大相位确定。   

                                  α与最大相位的关系:        sin(\phi _{m})=\frac{1-\alpha }{1+\alpha }

④现在找最大相位,超调量--阻尼比--相位裕度三者也有关系。举例超调量小于10%→阻尼比ζ就要大于0.6→相位裕度要大于60°。当前的相位裕度=10.4°,要提到60°,再来个50°就够了。再由于超前补偿器的特性,要额外多加几度,最终选定添加55°的相位超前。因此最大相位=55°。

⑤根据公式可计算得到α≈0.1

⑥根据下面关系式计算出超前补偿器所造成的幅度增加量

                                                          20lg(\frac{1}{\sqrt{\alpha }})\approx 10dB

⑦从(K=10的)bode图可得知,截止频率为3.49rad/s,幅值=-10dB的频率为6.1rad/s,因此引入α后新截止频率=6.1rad/s,利用这个新截止频率来计算T

w_{m}=\frac{1}{T\sqrt{\alpha }}\Rightarrow T=\frac{1}{w_{m}\sqrt{\alpha }}\approx 0.52

⑧至此超前补偿器的三个参数全部得到 K=10 α=0.1 T=0.52。看看波形图。

%重新将开环传递函数写一遍
s=tf('s');
P_pitch=(1.151*s+0.1774)/(s^3+0.739*s^2+0.921*s);

%引入超前补偿器
K=10;
alpha=0.1;
T=0.52;
C_lead=K*(T*s+1)/(alpha*T*s+1);
margin(C_lead*P_pitch),grid

%闭环输出响应
sys_cl=feedback(C_lead*P_pitch,1);
step(sys_cl),grid
title('加入了K α T后的拥有超前补偿器的输出阶跃响应图');

⑨输出波形图很难直观看到超调量等数据,采用stepinfo命令,可看到全部的性能指标值。

stepinfo(sys_cl)

⑩如果某些指标还是差一些,可以微调   K    α  的值,例如本例超调量依然略大,取α=0.0718(添加60°的相位超前)(T=0.5654),就可以解决问题。

滞后补偿器

与根轨迹里的滞后补偿器完全相同

适用条件:滞后补偿器可以在保持带宽频率不变的情况下增加低频段的幅值增益,这样就可以抬高稳定值,并且带宽频率不变

传递函数:                                             C(s)=\frac{s+z}{s+p}

零极点的选择:通常来说,与根轨迹法一样,滞后补偿器的零点要与系统的极点放在一起;至于极点选个靠近虚轴的就可以,例如 p=0.01 即可

%找到回路增益Kp,滞后补偿器C_lag后,命令如下(步骤同根轨迹)
bode(Kp*C_lag*P_pitch),grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%在控制系统设计器中也可以实现(步骤同根轨迹)
controlSystemDesigner('bode',P_pitch)

3.状态空间法

观察能控性(全状态反馈控制器)

先给出能控性矩阵M=[B  AB  A^2*B ....  A^(n-1)*B]

A = [-0.313 56.7 0; -0.0139 -0.426 0; 0 56.7 0];
B = [0.232; 0.0203; 0];
C = [0 0 1];
D = [0];

sys_order = order(motor_ss)   %找到最小的MIMO系统 

mo = ctrb(A,B);      %构建能控性矩阵M=[B AB ....]
sys_rank=rank(co)    %找到矩阵的秩


若sys_order=sys_rank,则矩阵满秩,则能控

观察能观性(状态观测器)

no=obsv(A,C);   %构建能观性矩阵N=[C CA ...]T
sys_rank=rank(no) 

设计全状态反馈控制器

经过对原状态空间进行全状态反馈的变换,可以得新状态空间表达式 

\dot{x}=(A-BK)x+Br

y=Cx

(1)做法一:极点放置法

适用于标准的一阶二阶系统(不适用于本例)

①从[sl-(A-BK)]的行列式(矩阵的根)中找到闭环极点,根据阶数可获得几个可任意放置的极点(一阶可放一个,二阶可放两个)

②先找到极点要放的位置,例如 p1=-5+i   p2=-5-i(举例),调用命令得到K值

%开环传递函数(举例)
J = 0.01;%转子的转动惯量
b = 0.1;%电机粘性摩擦常数   
K = 0.01;%电动势常数     电机转矩常数
R = 1;%电阻
L = 0.5;%电感


   %状态方程
A = [-b/J   K/J
      K/L   -R/L];
B = [0
	1/L];
C = [1   0];
D = 0;
motor_ss = ss(A,B,C,D)

   %放极点,求Kc
p1 = -5 + 1i;
p2 = -5 - 1i;
Kc = place(A,B,[p1 p2])

③Kc=   13.0100   -1.0000,将K带回新状态方程,看输出响应波形

motor_ss = ss(A,B,C,D)

motor_ss_new=ss(A-B*Kc,B,C,D)
t=0:0.01:3;
step(motor_ss_new,t)
title('状态控制器的输出响应')

④通常来说,放置极点不可能一次成功,一般都需要引入预补偿器Nbar,目的是缩放输入,抬高稳定值。

Nbar = rscale(motor_ss,Kc)%利用原状态方程算Nbar

t = 0:0.01:10;
step(motor_ss_new*Nbar,t)%将新状态方程×Nbar
grid
title('状态控制器+预补偿器的输出响应')

注意:rscale函数(查找比例因子以消除稳态误差函数)内容可以在matlab官方help处找到,复制粘贴到matlab当前目录中

(2)做法二:使用LQR(线性二次调节器)来生成最佳增益矩阵K (适用本例)

①若要使用LQR,需要定义两个参数:状态成本加权矩阵Q和控制权重矩阵R,为图简单,选择R=1,Q=p*C ' *C

②LQR参数全了,接下来得到反馈控制矩阵K,先令加权因子p=2

%状态空间方程
A = [-0.313 56.7 0;-0.0139 -0.426 0;0 56.7 0];
B = [0.232; 0.0203;0];
C = [0 0 1];
D = [0];
pitch_ss=ss(A,B,C,D);

%设计LQR
p=2;
Q=p*C'*C
R=1;
[K]=lqr(A,B,Q,R)

③可得K= -0.5034   52.8645    1.4142,根据此K得到新状态方程,且查看输出响应

%新状态方程
sys_cl=ss(A-B*K,B,C,D);

%查看输出响应
step(sys_cl),grid
title('LQR下的闭环阶跃响应')

④如果输出波形响应太慢,则增大LQR里的加权因子p值(本例将加权因子增大到p=50),这样可以得到更大的K,加快系统响应,但是与此同时也会降低稳态值。如果输出波形稳态值过小,稳态误差过大,则再调节LQR也没救了,需要再同样引入预补偿器Nbar

⑤同样使用rscale(查找比例因子以消除稳态误差函数)来得到Nbar

Nbar=rscale(A,B,C,D,K)%在LQR的基础上继续做,因此K为之前由LQR得到的K

⑥得到Nbar看看输出波形

%带LQR+Nbar的新状态空间方程
sys_cl=ss(A-B*K,B*Nbar,C,D);%也可以同做法一,先求得sys_cl_new,再×Nbar,是一样的

step(sys_cl),grid
title('LQR+Nbar下的闭环阶跃响应')

预补偿器Nbar

由于上述两种做法都可能用到Nbar,因此做一些补充

  •  添加的预补偿器N(-)并不位于反馈回路,这就意味着如果有未知干扰,预补偿器对此部分的偏差无法消除。
  • 换句话说,预补偿器的作用就是将输入放大或者缩小,导致了输出稳态值的提升,但是稳态误差说到底还是没真正被消除。根本原因就是N(-)不在反馈回路,它不能真正改变输出的波形,也就是说不会真正减小稳态误差
  • 有一种方法是将预补偿器与积分控制相结合,既可以缩放输出波形,又可以一定程度减小稳态误差

4.PID

PID传递函数:       

                                 C(s)=K_{p}+\frac{k_{i}}{s}+K_{d}s

%打开控制系统设计器
s=tf('s');
P_pitch=(1.151*s+0.1774)/(s^3+0.739*s^2+0.921*s);

controlSystemDesigner(P_pitch)

左上角--调节方法--PID调节--选择鲁棒响应时间/经典设计公式

  • 鲁棒响应时间就是自动谐调PID参数,可以为所有类型的PID调节参数,可用于设计稳定、不稳定、整合的系统。
  • 经典设计公式算法需要稳定的或带有积分的设备,无法调整导数滤波器,但是会有很多的公式方法。

此例选择鲁棒响应时间算法--设计模式时域表示--参数调调调--设计成功


5.数字控制器

(1)做法一:状态方程+LQR 

①先得到离散状态空间方程,使用c2d命令将连续时间模型→离散时间模型,c2d需要三个参数:系统模型+采样时间+保持电路的类型(通常用零阶保持器zoh)

②采样时间:采样频率至少比闭环bode图的带宽频率(闭环截止频率)(不是开环截止频率)大30倍。根据闭环波特图可得,带宽频率为2rad/sec,因此得到0.01s的采样时间。

③使用c2d命令,可以得到原离散状态空间模型

A = [-0.313 56.7 0;-0.0139 -0.426 0;0 56.7 0];
B = [0.232; 0.0203;0];
C = [0 0 1];
D = [0];
sys_ss=ss(A,B,C,D);
Ts=0.01;
sys_d=c2d(sys_ss,Ts,'zoh')

 ④验证能控性

co = ctrb(sys_d);
Controllability=rank(co)

⑤离散全状态反馈控制系统,得到新离散状态空间方程

x(k+1)=(A-BK)x(k)+Br(k)

y=Cx(k)

⑥依旧利用LQR,不过是离散版本

A=sys_d.A;
B=sys_d.B;
C=sys_d.C;
D=sys_d.D;
p=50;
Q=p*C'*C
R=1;
[K]=dlqr(A,B,Q,R)

⑦得到K,就可以得到离散状态空间方程了,且看看输出波形

time=0:0.01:10;
theta_des=ones(size(time));
sys_cl=ss(A-B*K,B,C,D,Ts);
[y,t]=lsim(sys_cl,theta_des);
%画图
stairs(t,y)
title('DLQR下的闭环阶跃响应')

⑧可能仍需要添加预补偿器Nbar,但是离散版本不可以用rescale自动求了,需要手动调调改改,最后选定Nbar=6.95。重新得离散状态空间方程+输出波形

Nbar=6.95;
sys_cl=ss(A-B*K,B*Nbar,C,D,Ts);%多乘了Nbar
[y,t]=lsim(sys_cl,theta_des);
%画图
stairs(t,y)
title('DLQR下的闭环阶跃响应')

(2)做法二:传递函数+根轨迹

①将连续传函转换成离散传函

Ts=1/50;
dp_pitch,=c2d(P_pitch,Ts,'zoh')
zpk(dp_pitch)%变成因式分解形式,更容易看零极点

②利用zgrid命令找到离散的 可接受符合要求的 根轨迹区域,然后在该区域里找增益K。而使用zgrid命令需要两个参数:自然频率+阻尼比。

w_{n}\geq \frac{1.8}{T_{r}}               \zeta =\frac{-ln(M_{p})}{\sqrt{\pi ^{2}+ln^{2}(M_{p})}}

假设上升时间为5秒,超调量为10%,因此固有频率大于0.36,阻尼比大于0.6。但由于zgrid命令的固有频率单位是rad/sample,因此换算成0.0072rad/sample

Wn=0.0072;
zeta=0.6;

%画根轨迹图
rlocus(dp_pitch)
zgrid(zeta,Wn)
axis([0.97 1 -0.05 0.05])

后面步骤与代码同根轨迹法 

③选极点→得Kp值→画闭环输出波形→不满意→乘上超前/滞后补偿器→重新选极点→得新Kp值→Kp×超前/滞后补偿器×原离散开环传递函数→得闭环传递函数→画闭环输出波形。

(3)做法三:传递函数+PID+根轨迹

①将连续传函转换成离散传函

Ts=1/50;
dp_pitch,=c2d(P_pitch,Ts,'zoh')
zpk(dp_pitch)

②把连续PID转换成离散PID

%离散PID
Kp = 100;
Ki = 200;
Kd = 10;
 
C = Kp + Ki/s + Kd*s;
dC = c2d(C,Ts,'tustin')

③用离散PID去控制离散传递函数,看看输出响应波形

sys_cl = feedback(dC*dP_pitch,1);
[x2,t] = step(sys_cl,12);  %12的意思是时间范围
%画图
stairs(t,x2)

如果发现输出响应不稳定,但该PID参数在连续域可以得到很好的调节结果(在离散域就调节不好),需要进根轨迹找原因

④看看离散系统+离散PID的根轨迹

%不稳定,看看根轨迹
rlocus(dC*dP_pitch)
axis([-1.5 1.5 -1 1])
title('离散系统带PID的根轨迹')

⑤经过很复杂的、我不会的分析后发现,PID控制器有个我们不想要的极点,因此就需要改变PID控制器.因此给离散PID增加了一个-0.82的极点(随便举例)

%给离散PID再增加一个极点
z = tf('z',Ts);
dC_new = dC/(z+0.82);
%看看新根轨迹
rlocus(dC_new*dP_motor)

⑥在新根轨迹上找个极点,看看对应的K值。此时得到K值

⑦在此K值下+新离散PID+原离散开环传递函数的输出波形如何?

sys_cl = feedback(K*dC_new*dP_pitch,1);%其中K值为上一步选择极点得到的值
[x3,t] = step(sys_cl,8);   %8的意思是时间范围

%画输出波形图
stairs(t,x3)
title('离散域下开传函+PID+根轨迹')

⑧如果不行就重新选K,一直到达到设计要求

四.附加资料

开环截止频率与闭环截止频率 - 知乎 (zhihu.com)

分贝dB与放大倍数的转换关系及对照表_db值对照表_李锐博恩的博客-CSDN博客

关于超前补偿器和滞后补偿器的理解_超前滞后补偿零极点位置的影响_zmhzmhzm的博客-CSDN博客

  • 9
    点赞
  • 79
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值