matlab 迭代pid参数自整定,求助：可以整定PID参数的MATLAB程序

我想采用可自整定PID参数的PID方法控制一个被控对象为：G(S)=2.45*e^(-8s)/[(10.5s+1)*(2.5s+1)]的传递函数，现在采用了一下RBF整定的程序，但是控制图像发散，不知道要怎么改，请高人指点。或者用其他的方法能够实现也可以，谢谢！

%Adaptive PID control based on RBF Identification

clear all;

close all;

ts=0.001;

sys=tf(0.093,[1,0.4952,0.038],'inputdelay',8);

dsys=c2d(sys,ts,'z');

[num,den]=tfdata(dsys,'v');

xite=0.25;

alfa=0.05;

belte=0.01;

x=[0,0,0]';

ci=30*ones(3,6);

bi=40*ones(6,1);

w=10*ones(6,1);

h=[0,0,0,0,0,0]';

ci_1=ci;ci_3=ci_1;ci_2=ci_1;

bi_1=bi;bi_2=bi_1;bi_3=bi_2;

w_1=w;w_2=w_1;w_3=w_1;

u_1=0;y_1=0;

u_2=0;y_2=0;

xc=[0,0,0]';

error_1=0;error_2=0;error=0;

%kp=rand(1);

%ki=rand(1);

%kd=rand(1);

kp=0.1;

ki=0.10;

kd=0.1;

kp_1=kp;

kd_1=kd;

ki_1=ki;

xitekp=0.5;

xitekd=0.5;

xiteki=0.5;

for k=1:1:2000

time(k)=k*ts;

rin(k)=1;

yout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2;  %Nonlinear plant

for j=1:1:6

h(j)=exp(-norm(x-ci_1(:,j))^2/(2*bi_1(j)*bi_1(j)));

end

ymout(k)=w_1'*h;

d_w=0*w;      % Defining matrix number of d_w equal to that of w

for j=1:1:6

d_w(j)=xite*(yout(k)-ymout(k))*h(j);

end

w=w_1+d_w+alfa*(w_1-w_2)+belte*(w_2-w_3);

d_bi=0*bi;

for j=1:1:6

d_bi(j)=xite*(yout(k)-ymout(k))*w_1(j)*h(j)*(bi_1(j)^-3)*norm(x-ci_1(:,j))^2;

end

bi=bi_1+ d_bi+alfa*(bi_1-bi_2)+belte*(bi_2-bi_3);

for j=1:1:6

for i=1:1:3

d_ci(i,j)=xite*(yout(k)-ymout(k))*w_1(j)*h(j)*(x(i)-ci_1(i,j))*(bi_1(j)^-2);

end

end

ci=ci_1+d_ci+alfa*(ci_1-ci_2)+belte*(ci_2-ci_3);

%%%%%%%%%%%%%%%%%%%%%%Jacobian%%%%%%%%%%%%%%%%%%%%%%%

yu=0;

for  j=1:1:6

yu=yu+w(j)*h(j)*(-x(1)+ci(1,j))/bi(j)^2;

end

dyout(k)=yu;

%%%%%%%%%%%%%%%%%%%%%%Start of Control system%%%%%%%%%%%%%%%%%%

error(k)=rin(k)-yout(k);

kp(k)=kp_1+xitekp*error(k)*dyout(k)*xc(1);

kd(k)=kd_1+xitekd*error(k)*dyout(k)*xc(2);

ki(k)=ki_1+xiteki*error(k)*dyout(k)*xc(3);

if kp(k)<0

kp(k)=0;

end

if kd(k)<0

kd(k)=0;

end

if ki(k)<0

ki(k)=0;

end

M=1;

if M==2  %Only PID Control

kp(k)=0.03;

ki(k)=0.01;

kd(k)=0.03;

end

u(k)=u_1+kp(k)*xc(1)+kd(k)*xc(2)+ki(k)*xc(3);

if u(k)>=10

u(k)=10;

end

if u(k)<=-10

u(k)=-10;

end

%Return of parameters

x(1)=u(k);

x(2)=yout(k);

x(3)=y_1;

u_1=u(k);

u_2=u_1;

y_1=yout(k);

y_2=y_1;

w_3=w_2;

w_2=w_1;

w_1=w;

xc(1)=error(k)-error_1;             %Calculating P

xc(2)=error(k)-2*error_1+error_2;   %Calculating D

xc(3)=error(k);                     %Calculating I

error_2=error_1;

error_1=error(k);

kp_1=kp(k);

kd_1=kd(k);

ki_1=ki(k);

end

figure(1);

plot(time,rin,'b',time,yout,'r');

xlabel('time(s)');ylabel('rin,yout');

figure(2);

plot(time,yout,'r',time,ymout,'b');

xlabel('time(s)');ylabel('yout,ymout');

figure(3);

plot(time,dyout);

xlabel('time(s)');ylabel('Jacobian value');

figure(4);

subplot(311);

plot(time,kp,'r');

xlabel('time(s)');ylabel('kp');

subplot(312);

plot(time,ki,'r');

xlabel('time(s)');ylabel('ki');

subplot(313);

plot(time,kd,'r');

xlabel('time(s)');ylabel('kd');

[本帖最后由 leoyfli 于 2009-6-7 22:21 编辑]