Octave实现位置式PID算法

由于Matlab不让用，只能“你不让爷用，爷就用别的”，选择开源的Octave以及scilab进行相关领域的学习。Octave的代码和Matlab几乎是100%相同的，只有一些专用的包的函数，可能有些还没来得及写，或者有些差异。但这种差异，新手一般体会不到，老手应该能自己解决了吧。

数字PID控制

PID控制器是一种线性控制器，根据给定值 $r_i$ 与实际值构成控制偏差 $err(t)=r_i(t)-y_o(t)$ 。其控制规律为

$u(t)=k_p(err(t)+\frac{1}{T_1}\int^{t}_{0}err(t)\text dt+\frac{T_D\text derr(t)}{\text dt})$

写成传递函数的形式为

$G(s)=\frac{U(s)}{E(s)}=k_p(1+\frac{1}{T_1s}+T_Ds)$

其中， $k_p$ 为比例系数； $T_1$ 为积分时间常数； $T_D$ 为微分时间常数。 $U (s)$ 为输出量的拉氏量， $E (s)$ 为输入量的拉氏量。

PID控制器各校正环节如下：

比例环节：成比例地反应控制系统的偏差信号 $err (t)$ ，偏差一旦产生，控制器立即产生控制作用，从而减少偏差。
积分环节：用于消除静差，提高系统的无差度。积分作用的强弱与积分时间常数 $T_1$ 成负相关。
微分环节：反应偏差信号的变化趋势（变化速率），并能在偏差信号变得太大之前，在系统中引入一个有效的早期修正信号，从而加快系统的动作速度，减少调节时间。

位置式PID控制算法

以一系列的采样时刻点 $k T$ 代表连续时间 $t$ ，以举行法数值积分代替积分，以一阶后向差分近似代替微分，即

$\left\{\begin{aligned} &t\approx kT, k=0,1,2...\\ &\int^1_0err(t)\text dt\approx T\sum^k_{j=0}err(jT)=T\sum^k_{j=0}err_j\\ &\frac{\text derr(t)}{\text dt}\approx\frac{err(kT)-err((k-1)T)}{T}=\frac{err_k-err_{k-1}}{T} \end{aligned}\right.$

可得离散PID表达式

$\begin{aligned} u(k)=&k_p(err_k+\frac{T}{T_1}\sum^k_{j=0}err_j+\frac{T_D}{T}(err_k-err_{k-1}))\\ =&k_perr_k+k_I\sum^k_{j==0}err_jT+k_d\frac{err_k-err_{k-1}}{T} \end{aligned}$

其中， $k_I=\frac{k_p}{T_1},k_d=k_pT_D$ ， $T$ 为采样周期， $k$ 为采样序号， $k = 1, 2, ...$ 。其控制系统为

在这里插入图片描述

偏微分方程求解

设被控对象为电机模型传递函数

$G(s)=\frac{1}{Js^2+Bs},J=0.0067,B=0.10$

则我们可以通过Octave通过ODE45的方法来求解方程，输入信号为 $r_{in}(k)=0.50\sin(2\pi t)$ ，采用PID控制方法设计控制器，其中 $k_p=20.0,k_d=0.50$ 。

其中，ODE45的调用方法为

[t,x]=ode45(func,tspan,x0,op,para)

其返回值t是一个列向量，x是一个矩阵；参数func为待处理函数或其路径，tspan=[t0 tf]为微分方程组的积分区间，

代码为

ts = 0.001; %Sampling time
xk = zeros(2,1);
e1 = 0; u1 = 0; %初始化误差

time = ts:ts:2000*ts
rin = 0.5*sin(1*2*pi*time);
for k = 1:1:2000
  para = u1;
  tSpan = [0 ts];
  [tt, xx] = ode45("test1_func", tSpan, xk, [], para);
  xk = xx(length(xx),:);
  yout(k) = xk(1);
  e(k) = rin(k) - yout(k);  %误差
  de(k) = (e(k)-e1)/ts;     %误差的一阶导数
  u(k) = 20.0*e(k)+0.5*de(k);
  u(k) = min(u(k),10.0);    %限制u(k)所在区间为[-10,10]
  u(k) = max(u(k),-10.0);
  u1 = u(k);
  e1 = e(k);
end

subplot(1,2,1)
plot(time, rin, 'r', time ,yout, 'b');
xlabel('time(s)'), ylabel('rin,yout');

subplot(1,2,2)
plot(time, rin-yout, 'r');
xlabel('time(s)'), ylabel('error');

函数文件为

% test1_func.m
function dy=PlantModel(t,y,flag,para)
  u = para;J = 0.0067;B = 0.1;
  
  dy = zeros(2,1);
  dy(1) = y(2);
  dy(2) = -(B/J)*y(2)+(1/J)*u;

得到其结果为

在这里插入图片描述
此外，可以通过XCOS进行仿真，被控对象为三阶传递函数，采用XCOS与脚本结合的方式。主程序由XCOS实现，控制由scilab实现。

离散系统的数字PID控制仿真

控制对象为

$G(S)=\frac{523500}{s^3+87.35s^2+10470s}$

采样时间为1ms，采用z变换进行离散化，经过z变换后的离散化对象为

$\begin{aligned} y_{out}=&-den(2)y_{out}(k-1)-den(3)y_{out}(k-2)-den(4)y_{out}(k-3)\\ &+num(2)u(k-1)+num(3)u(k-2)+num(1)u(k-3) \end{aligned}$

其Octave控制代码为

ts = 0.001; %采样时间
N = 1000;    %采样个数

pkg load control    %载入control包，内含tf函数
pkg load tisean     %内含c2d函数

sys = tf(2.235e5,[1,87.35,1.047e4,0]); 
dsys = c2d(sys,ts,'z'); %z变换
[num,den] = tfdata(dsys, 'v');

titles = ["step signal";"square signal";"sin signal"]

for s = 1:3
%参数初始化
u = zeros(1,3);y = zeros(1,3);x = zeros(1,3);
err = 0;

time = ts:ts:ts*N;
rin = ones(1,N);
PDI = [2 0.001 0.001];          %表示P，D，I的分量
if s==2
  rin = sign(sin(2*2*pi*time)); %方波信号
elseif s==3
  rin = 0.5*sin(2*2*pi*time);   %正弦信号
endif

for k = 1:1:N
  yout(k) = -sum(den(2:4).*y)+sum(num.*u);
  error(k) = rin(k)-yout(k);

  u(2:3) = u(1:2);  
  u(1) = sum(PDI.*x);  %PID 控制
  u(1) = sort([u(1),10,-10])(2)   %区间限制
  
  y = [yout(k) y(1:2)];
  
  x(1) = error(k);
  x(2) = (error(k)-err)/ts;
  x(3) = x(3) + error(k)*ts;
  
  err = error(k);
end

subplot(230+s)
plot(time, rin, 'b', time, yout, 'r')
xlabel('time(s)'), ylabel('rin, rout')
title(titles(s,:))

subplot(233+s)
plot(time,error)
xlabel('time(s)'),ylabel('error')
end