程式碼如下:
clear
Ts=25;
p=1.00151*10^(-4);
q=8.67973*10^(-3);
r=40.0;
a=exp(-p*Ts);
b=(q/p)*(1-exp(-p*Ts));
yo =25;
y(1)=25;
for k=1:500
if (k<=40)
ref(k)=35;
elseif ((k>40)& (k<=80))
ref(k)=50;
elseif ((k>80)& (k<=120))
ref(k)=65;
elseif (k>120)
ref(k)=80;
end
GE=1/15; GC=1/15;
GU=450;
for k=1:1:500
e(k)=ref(k)-y(k);
ee(k)=e(k)*GE;
if (k==1)
ce(k)=0;
else
ce(k)=e(k)-e(k-1);
end
cee(k)=e(k)*GC;
mu(k)=1;
u(k)=mu(k)*GU;
y(k+1)=a*y(k)+b/(1+exp(0.5*y(k)-r))*u(k)+(1-a)*yo
end
subplot(2,2,1)
plot(0:25:500* Ts,y(1), 0:25:500* Ts, ref,’