mohuoidobj.asv

氧乐果控制过程

function obj=mohupidobj(chrom)
global rin yout timef tu
ts=20;

%Delay plant
sys=tf([1],[60,1],'inputdelay',80);
dsys=c2d(sys,ts,'zoh');
[num,den]=tfdata(dsys,'v');
[n1,n2]=size(chrom)
n=0;
a1=0;
while (n<n1)
n=n+1;
error_1=0;error_2=0;
ei=0;
rin=1.0;
u_1=0;u_2=0;u_3=0;u_4=0;u_5=0;
y_1=0;y_2=0;y_3=0;
x=[0,0,0]';
B=0;
error_1=0;
tu=1;
s=0;
P=100;

%------------ Return of PID parameters -------------
for k=1:1:P
r(k)=rin;
timef(k)=k*ts;
yout(k)=-den(2)*y_1+num(2)*u_5;
error(k)=r(k)-yout(k);
u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k);   
y_3=y_2;y_2=y_1;y_1=yout(k);
x(1)=abs(error(k)/rin);              %compute e
x(2)=abs((error(k)-error_1)/(ts*rin)); %compute ec
error_2=error_1;
error_1=error(k);

%对e,ec进行量化
 if x(1)>=0.8 
e=5*sgn(error(k))
elseif x(1)>=0.5
e=4*sgn(error(k))
elseif x(1)>=0.3
e=3*sgn(error(k))
elseif x(1)>=0.1
e=2*sgn(error(k))
elseif x(1)>=0.03
e=sgn(error(k))
elseif
e=0
 end

 if x(2)>=0.8 
ec=5*sgn(error(k))
 elseif x(2)>=0.5
ec=4*sgn(error(k))
 elseif x(2)>=0.3
ec=3*sgn(error(k))
 elseif x(2)>=0.1
ec=2*sgn(error(k))
elseif x(2)>=0.03
 ec=sgn(error(k))
elseif
 ec=0
  end

u(k)=a1*e+(1-a1)*ec; 
%    if u(k)>=10
%       u(k)=10;
%    end
%    if u(k)<=-10
%       u(k)=-10;
%    end   
%    

   if s==0
%    if yout(k)>0.95&yout(k)<1.05
%       tu=timef(k);
%        s=1;
%    end 
     if yout(k)>0.99
       tu=timef(k);
        s=1;
     end
 end
end

for i=1:1:P
% Ji(i)=0.999*abs(error(i))+0.05*u(i)^2;
% Ji(i)=abs(error(i));
% Ji(i)=(abs(error(i))*ts)^2;
% Ji(i)=(abs(error(i))*ts);
   B=B+Ji(i);   
  if i>1   
   erry(i)=yout(i)-yout(i-1);
   if erry(i)<0
  B=B+59*abs(erry(i));
%    B=B+10*abs(erry(i))
   end    
  end
end

% obj(n)=B+0.2*tu*10;
 obj(n)=B;
end
obj=obj';
 
% plot(obj)