%% AFC控制器
%溶解氧利用AFC控制,无监督控制器,无参数投影
%氧气转换系数控制范围是0-240/d,内回流量范围是0-92230.即5*Qin0
%% 清空
clc;
clear all;
close all;
%% BSM1初始化
openloopinit;
step_day = 0.0125;
timestep = step_day/24;%0.01;
timeend =1; %30; %步长0.01,仿真30天 效果可以
xo1 = XINIT1(1,1:13); %第一分区
xo2 = XINIT2(1,1:13); %第二分区
xo3 = XINIT3(1,1:13); %第三分区
xo4 = XINIT4(1,1:13); %第四分区
xo5 = XINIT5(1,1:13); %第五分区
xoSettler = SETTLERINIT; %二沉池
Z1His = []; %第一分区的组分历史
Z2His = []; %第二分区的组分历史
Z3His = []; %第三分区的组分历史
Z4His = []; %第四分区的组分历史
Z5His = []; %第五分区的组分历史
Z_SettlelHis = [];%二沉池组分历史
Z_CompEffluentHis = [];%二沉池出水组分历史
BOD_5e_His = []; %出水BOD历史
COD_e_His = []; %出水COD历史
TSS_Settler = 6394; %TSS_Settler为二沉池固体悬浮物浓度
Qr = Qin0; %Qr是污泥回流量,从二沉池回流到第一分区,Qin0=18446
Qw = 385; %Qw是剩余污泥排放量
%Unit5ToUnit1,污水内回流组分,从第五分区流到第一分区,13种组分,TSS5为第五分区固体悬浮物浓度,Qintr为内回流量=3*Qin0
Unit5ToUnit1 = [ S_I5 S_S5 X_I5 X_S5 X_BH5 X_BA5 X_P5 S_O5 S_NO5 S_NH5 S_ND5 X_ND5 S_ALK5 TSS5 Qintr];
%SettlerToUnit1,二沉池污泥回流组分,从二沉池到生化池第一分区,TSS_Settler为二沉池固体悬浮物浓度
SettlerToUnit1 = [30,0.88949,2247.0367,96.4143,5004.6489,292.9183,884.2618,0.49094,10.4152,1.7334,0.68828,6.8972,4.1256,TSS_Settler, Qr];
CompEffluent = zeros(1,15);%出水组分
%% 溶解氧PID控制器参数
Exp_DoUnit5 = 2; %溶解氧浓度设定值为2mg/L
Exp_DoUnit5M=[]; %记录设定值线
Do5_His=[]; %记录BSM1模型真实DO值
Err_Do5_His = []; %记录第五分区溶解氧误差
KLa5_His = []; %记录第五分区的曝气量
Error_Do5 = 0; %溶解氧的误差
preErr_Do5 = 0;
pre2Err_Do5 = 0;
Kp_Do = 200;
Ti_Do = 15;
Td_Do = 2;
% Kp_Do = 10;
% Ti_Do = 2;
% Td_Do = 0.5;
%% 硝态氮PID控制器参数
Exp_SNoUnit2 = 1; %硝态氮浓度设定值为1mg/L
Exp_SNoUnit2M=[]; %用于画图,画设定值线
SNo2_His=[]; %记录BSM1模型真实硝态氮值
Err_SNo2_His = []; %记录第二分区硝态氮误差
Qintr_His = []; %记录内回流量
Error_SNo2 = 0;
preErr_SNo2 = 0;
pre2Err_SNo2 = 0;
Kp_SNo2 = 50000; %三个参数分别为10000 ; 3000 ; 100 时效果不错
Ti_SNo2 = 5000;
Td_SNo2 = 400;
%% 上层指标
AE = 0; %曝气能耗
PE = 0; %泵送能耗
EffQuality = 0; %出水水质
%% 执行控制
theta_0=[10 10 10]; %初始控制率!!!任意设定
abs_theta_his=[]; %用于记录abs_theta的值,看他是否有界
to = clock; %计算运行时间
epoch = 1; %每一个控制周期,epoch+1
for time = 0:timestep:timeend-2*timestep
epoch
% %第五分区溶解氧阶跃变化
if rem(epoch,2/step_day) == 0
Exp_DoUnit5 = 0.2+(2-0.2)*rand(1); %0.4-3 溶解氧优化值变化范围
Exp_SNoUnit2 = 0.2+(2-0.2)*rand(1); %0.2-2 硝态氮优化值变化范围
end
%% 载入BSM1数据
% Flow_comb1 = combiner([CONSTINFLUENT(floor(epoch/(0.25/step_day))+1,2:16) SettlerToUnit1]);%输入组分和二沉池中回流的组分(1*30)混合以后的输出 1*15
Flow_comb1 = combiner([DRYINFLUENT(floor(epoch/(0.25/step_day))+1,2:16) SettlerToUnit1]);%输入组分和二沉池中回流的组分(1*30)混合以后的输出 1*15
% Flow_comb1 = combiner([RAININFLUENT(floor(epoch/(0.25/step_day))+1,2:16) SettlerToUnit1]);%输入组分和二沉池中回流的组分(1*30)混合以后的输出 1*15
% Flow_comb1 = combiner([STORMINFLUENT(floor(epoch/(0.25/step_day))+1,2:16) SettlerToUnit1]);%输入组分和二沉池中回流的组分(1*30)混合以后的输出 1*15
Flow_comb2 = combiner([Flow_comb1 Unit5ToUnit1]); %%混合以后与内回流组分再次混合1*30后的输出 1*15
unit1_In = carboncombiner(CARBONSOURCECONC,[carb1 Flow_comb2]);%第一分区预输入,与外加碳源混合 1*15
CompInput1 = [unit1_In KLa1]; %第一分区的输入,1*16
dx1dt = ASM1( xo1 ,CompInput1 ,VOL1);
xo1 = xo1 + dx1dt'*timestep; %
for i = 1:13
if xo1(i) < 0.0
xo1(i) = 0.0;
else
xo1(i) = xo1(i);
end
end
Xf1=X_I2TSS*xo1(3)+X_S2TSS*xo1(4)+X_BH2TSS*xo1(5)+X_BA2TSS*xo1(6)+X_P2TSS*xo1(7);
unit1_Out = [xo1 Xf1 CompInput1(15)]; %第一分区的输出组分浓度,Xf为Xi,TSS;Q1为流量
unit2_In = carboncombiner(CARBONSOURCECONC,[carb2 unit1_Out]);
CompInput2 = [unit2_In KLa2];
dx2dt = ASM1( xo2 ,CompInput2 ,VOL2);
xo2 = xo2 + dx2dt'*timestep; %
for i = 1:13
if xo2(i) < 0.0
xo2(i) = 0.0;
else
xo2(i) = xo2(i);
end
end
Xf2=X_I2TSS*xo2(3)+X_S2TSS*xo2(4)+X_BH2TSS*xo2(5)+X_BA2TSS*xo2(6)+X_P2TSS*xo2(7);
unit2_Out = [xo2 Xf2 CompInput2(15)];%第二分区的输出组分浓度
unit3_In = carboncombiner(CARBONSOURCECONC,[carb3 unit2_Out]);
CompInput3 = [unit3_In KLa3];
dx3dt = ASM1( xo3 ,CompInput3 ,VOL3);
xo3 = xo3 + dx3dt'*timestep; %
for i = 1:13
if xo3(i) < 0.0
xo3(i) = 0.0;
else
xo3(i) = xo3(i);
end
end
Xf3=X_I2TSS*xo3(3)+X_S2TSS*xo3(4)+X_BH2TSS*xo3(5)+X_BA2TSS*xo3(6)+X_P2TSS*xo3(7);
unit3_Out = [xo3 Xf3 CompInput3(15)];%第三分区的输出组分浓度
unit4_In = carboncombiner(CARBONSOURCECONC,[carb4 unit3_Out]);
CompInput4 = [unit4_In KLa4];
dx4dt = ASM1( xo4 ,CompInput4 ,VOL4);
xo4 = xo4 + dx4dt'*timestep; %
for i = 1:13
if xo4(i) < 0.0
xo4(i) = 0.0;
else
xo4(i) = xo4(i);
end
end
Xf4=X_I2TSS*xo4(3)+X_S2TSS*xo4(4)+X_BH2TSS*xo4(5)+X_BA2TSS*xo4(6)+X_P2TSS*xo4(7);
unit4_Out = [xo4 Xf4 CompInput4(15)];%第四分区的输出组分浓度
unit5_In = carboncombiner(CARBONSOURCECONC,[carb5 unit4_Out]);
CompInput5 = [unit5_In KLa5];
dx5dt = ASM1( xo5 ,CompInput5 ,VOL5);
xo5 = xo5 + dx5dt'*timestep; %
for i = 1:13
if xo5(i) < 0.0
xo5(i) = 0.0;
else
xo5(i) = xo5(i);
end
end
Xf5=X_I2TSS*xo5(3)+X_S2TSS*xo5(4)+X_BH2TSS*xo5(5)+X_BA2TSS*xo5(6)+X_P2TSS*xo5(7);
unit5_Out = [xo5 Xf5 CompInput5(15)];%第五分区的输出组分浓度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------溶解氧控制器作用-------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%第一步:参数初始化
Error_Do5 = Exp_DoUnit5 - xo5(8);% 当前溶解氧误差,溶解氧设定值-BSM1模型输出溶解氧真实值
KLa5_max=240; %控制量的最大值
k=1; %k为任意正定矩阵
gama=200; %自适应参数
E=[Error_Do5]'; %模糊系统输入量,即溶解氧误差
Q=[10]; %Q为任意n*n正定矩阵
A=[-k];
P=lyap(A',Q); %利用Lyapunov方法求出P=5
%第二步:主控制器包含3条模糊规则
%第1条模糊规则
a1=0.05;b1=-0.1;
if Error_Do5>=a1
uc(1)=0;
end
if Error_Do5>=(a1+b1)/2 && Error_Do5<=a1
uc(1)=2*( (Error_Do5-a1)/(b1-a1) )^2;
end
if Error_Do5>=b1 && Error_Do5<=(a1+b1)/2
uc(1)=1-2*( (Error_Do5-b1)/(b1-a1) )^2;
end
if Error_Do5<=b1
uc(1)=1;
end
%第2条模糊规则
uc(2)=exp ( -(Error_Do5 -0.00)^2 / (2*0.02^2) ); %第2条隶属函数
%第3条模糊规则
a3=-0.05; b3=0.1;
if Error_Do5<=a3
uc(3)=0;
end
if Error_Do5<=(a3+b3)/2 && Error_Do5>=a3
uc(3)=2*( (Error_Do5-a3)/(b3-a3) )^2;
end
if Error_Do5<=b3 && Error_Do5>=(a3+b3)/2
uc(3)=1-2*( (Error_Do5-b3)/(b3-a3) )^2;
end
if Error_Do5>=b3
uc(3)=1;
end
%第三步:模糊系统输出
FS1=0; %初始化为0,用于下面的累加
for i=1:1:3
FS2(i)=uc(i); %FS2为模糊系统的分子
FS1=FS1+uc(i); %FS1为模糊系统的分母
end
FS=FS2/FS1; %FS模糊系统的输出
KLa5_Uc=theta_0*FS' %计算主控制器输出
%第四步:参数投影
SS=E'*P*theta_0*FS';
abs_theta=norm(theta_0,2); %求theta的2范数
if abs_theta<KLa5_max||(abs_theta==KLa5_max && SS<=0)
dtheta=gama*E'*P*FS %自适应律保持不变
elseif abs_theta>= KLa5_max && SS>0
dthe
a_demo_AFC1_jar7ry_matlab_DEMO_matlabafc函数_AFC_源码
版权申诉
5星 · 超过95%的资源 8 浏览量
2021-10-04
02:40:40
上传
评论
收藏 6KB ZIP 举报
耿云鹏
- 粉丝: 62
- 资源: 4760
评论1