神经模糊技术的智能过程控制
本篇论文中采用由等效的四层链接网络构成的神经模糊控制器(NFC)作为过程反馈的控制器，利用学习算法，NFC能够学习通过更新模糊规则和隶属函数来自适应地控制过程。本篇用NFC及MNN实现对开环不稳定非线性CSTR的直接自适应控制。
- NFC神经网络：
分析：
NFC神经网络分为四层；
第一层：输入层，X1是变换过程输出误差；X2是变换后的误差变换。
第二层：语言层，用高斯函数作为隶属函数来确定观察到的信号。
第三层：规则层，实现与先决条件相关联的连接，其规则是每个规则节点只有一个来自语言变量的线性连接。
第四层：输出层，所有结果连接都完全连接到输出节点，视为输出强度。
NFC神经网络各层输入输出函数
程序实现如下：
function [ yt ] = NFC (u1,u2,a,b,w) %定义NFC函数，输入变量是u，输出变量是yt
global O1 O2 O3 k h; % 01代表第一层，02第二层，03第三层
%设定NFC结构为2-14-49-1，即2个输入，14个语言变量，49条控制规则，1个输出
%搭构神经模糊控制器的结构
%第一层 输入层为两个输入
u = [u1,u2];
for i=1:2
for j=1:7
O1(i,j) = u(i);
end
end
%第二层
for i=1:2
for j=1:7
I2 = - (O1(i,j) - a(i,j))^2/(b(i,j)^2);
O2(i,j) = exp(I2);
end
end
%第三层
for i=1:7
for j=1:7
I3 = O2(1,i)*O2(2,j);
O3(i,j) = I3;
end
end
I4 = 0; %定义
sum = 0;
%第四步
for i=1:1:7
for j=1:1:7
I4 = I4 + O3(i,j)*w(i,j);
sum = sum + O3(i,j);
end
end
yt = I4/sum;
end
- MNN估计器：
MNN估计器采用三层算法，第一层输入层，第二层隐藏层，第三层输出层，每层定义如下：
3.开环不稳定非线性CSTR实现：
- 求解该例中的初值，其微分方程如下:微分方程用“四阶龙格库塔”搭建模型求解，其程序如下：
function dotstate=runge_kutta(~,state)
global u k h;
t=k*h;
dotstate(1:2,1) = zeros(2,1);
Da = 0.072;
phi = 20;
B = 8;
delta = 0.3;
dotstate(1) = - state(1) + Da*(1-state(1))*exp(state(2)/(1+state(2)/phi));
dotstate(2) = -(1+delta)*state(2) + B*Da*(1-state(1))*exp(state(2)/(1+state(2)/phi))+delta*u(k);
end
- 实现主程序：
clear
clf
global u O1 O2 O3 k h; %定义全局变量
%初始 隶属度函数分布情况说明
a(1:2,:) = [ linspace(-1,1,7)
linspace(-1,1,7)]; %高斯函数的中心，从-1到1平均分配
b(1:2,1:7) = 0.25; %高斯函数的宽度参数
%权重赋予初值
w = [ -1 -1 -2/3 -2/3 -1/3 -1/3 0
-1 -2/3 -2/3 -1/3 -1/3 0 1/3
-2/3 -2/3 -1/3 -1/3 0 1/3 1/3
-2/3 -1/3 -1/3 0 1/3 1/3 2/3
-1/3 -1/3 0 1/3 1/3 2/3 2/3
-1/3 0 1/3 1/3 2/3 2/3 1
0 1/3 1/3 2/3 2/3 1 1 ];
%初值 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------一次循环结束后 将上次的结果进行转换储存-------
m1 = 4;
m2 = 5;
etax = 0.3;
bx = 0.01;
ax = 1;
w2x = zeros(5,4);
w2x(:,:) = 0.02*(rand(m2,m1)-0.5);
theta2x = 0.01*(rand(m2,1)-0.5);
w3x = zeros(5);
w3x = 0.02*(rand(m2,1)-0.5);
theta3x = 0.02*(rand(1)-0.5);
ax1 = ax;
w21 = w2x;
theta2x1 = theta2x;
w31 = w3x;
theta3x1 = theta3x;
ax0 = ax1;
w20 = w21;
theta2x0 = theta2x1;
w30 = w31;
theta3x0 = theta3x1;
s1(1:m1) = 0;
s2(1:2) = 0;
yx = 0;
q = 2;
eta = 0.8;
beita = 0.01;
y(1,1) =0;
u(1,1) =0;
%可以不定义变量，避免循环中给扩充矩阵而定义
delta2(1:5) = 0;
dE_dw(1:7,1:7) = 0;
dE_da(1:2,1:7) = 0;
dE_db(1:2,1:7) = 0;
x01 = 0.144;
x02 = 0.886; %平衡点A的初值
x(1:2,1) = [x01 x02]; % 模型（微分方程）的初值
total =40; %仿真时间
h =0.01; %步长
yd =2.750; %期望值
e(1:2,1) =[yd 0]; % e 和 ce 的初
k3 = 5; %设定k3的初值
%系统循环迭代推理
for k=1:1:(total/h) %(k=1:h:total)
% 对 e 和 ce 进行非线性饱和约束，将误差和误差变化控制在规则表之内，并量化(模糊化)映射到论域[-1.1]内部
e1 = (1 - exp(-10*e(1,k)))/(1 + exp(-10*e(1,k)));
e2 = (1 - exp(-0.01*e(2,k)))/(1 + exp(-0.01*e(2,k)));
%由NFC 得到系统的控制率
u(k) = k3 * NFC(e1,e2,a,b,w);
%求微分方程,四阶龙格-库塔法求解微分方程（系统模型） x为系统的状态量。
c1=runge_kutta(k*h,x(:,k));
c2=runge_kutta(k*h+h/2,x(:,k)+h*c1/2);
c3=runge_kutta(k*h+h/2,x(:,k)+h*c2/2);
c4=runge_kutta(k*h+h,x(:,k)+h*c3);
x(:,k+1)=x(:,k)+h*(c1+2*c2+2*c3+c4)/6;
%记录输出的误差
e(1,k+1) = yd - x(2,k);
e(2,k+1) = e(1,k+1) - e(1,k);
%系统输出()
y(k) = x(2,k);
%%%%------------------MNN 各层网络值的计算----------------------%%%%
if k<3
s1(1:4) = [y(1),0,u(1),0];
else
for j=1:1:m1 %MNN step1（m1=4）
if j < 3
s1(j) = y(k-j+1);
else
s1(j) = u(k-j+1+q); %q=2
end
end
end
sum2 = 0;
sum3 = 0;
for i = 1:1:m2 %MNN step 2 (m2=5)
for j = 1:1:m1
sum2 = sum2+w2x(i,j)*s1(j);
end
net2 = sum2 + theta2x(i);
s2(i) = (1-exp(-net2))/(1+exp(-net2));
sum3 = sum3 +w3x(i)*s2(i); % MNN step 3
end
net3 = sum3 + theta3x;
yx = ax*(1-exp(-net3))/(1+exp(-net3)); % MNN 结果，求出y的估计值。
%%%------------MNN中 求y对u* 的导数值------------------%%%
% y对u* 的导数
delta3 = (1/2)*ax*(1-yx/ax)*(1+yx/ax);
for i=1:1:m2
delta2(i) = 0.5*(1-s2(i))*(1+s2(i));
end
%-----
dy_du = 0;
for i=1:1:m2
dy_du = dy_du + w3x(i)*delta2(i)*w2x(i,q+1);
end
dy_du = dy_du*k3*delta3;
%%%-----------MNN 求取结果以后 利用梯度下降法更新 各层网络权值 和 阈值 以及特殊值-------------%%%
%%%%%%%%%%---------- 保存 k 时刻的值-----------------------
ax1 = ax;
w21 = w2x;
theta2x1 = theta2x;
w31 = w3x;
theta3x1 = theta3x;
%%%%%%%%%%%%%%%%%%%% %更新 w2 w3 theta2 theta3 ax %%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:1:m2
for j=1:1:m1
w2x(i,j) = w2x(i,j)+etax*(y(k)-yx)*delta3*w3x(i)*delta2(i)*s1(j)+bx*(w2x(i,j)-w20(i,j));
end
theta2x(i) = theta2x(i) + etax*(y(k)-yx)*delta3*w3x(i)*delta2(i)+bx*(theta2x(i)- theta2x0(i));
w3x(i) = w3x(i,1)+etax*(y(k)-yx)*delta3*s2(i)+bx*(w3x(i)-w30(i));
end
theta3x = theta3x + etax*(y(k)-yx)*delta3 +bx *(theta3x - theta3x0);
ax = ax + etax*(y(k)-yx)*(yx/ax) + bx*(ax-ax0);
%%%%%%%%% %将k时刻保存的值赋给下一个循环中作为：（k-1）时刻的值
ax0 = ax1;
w20 = w21;
theta2x0 = theta2x1;
w30 = w31;
theta3x0 = theta3x1;
% 初值 ********************************************************
% eta = 0.8;
% beita = 0.01；
%%% 求E对w的偏导---------------------------------
O3_sum = sum(sum(O3));
for i=1:1:7
for j=1:1:7
dE_dw(i,j) = -(yd -y(k))*dy_du*O3(i,j)/O3_sum;
end
end
O2_sum = sum(O2);
O3w_sum = 0;
for i=1:1:7
for j=1:1:7
O3w_sum = O3(i,j)*w(i,j);
end
end
SO2a1 = 0;
SO2a2 = 0;
SO2b1 = 0;
SO2b2 = 0;
for l = 1:1:7
SO2a1 = SO2a1 + O2(2,l);
SO2a2 = SO2a2 + O2(1,l);
SO2b1 = SO2b1 + O2(2,l);
SO2b2 = SO2b2 + O2(1,l);
end
%计算 dE/da1 dE/da2 dE/db1 dE/db2;
dEa1_sum = 0;
dEa2_sum = 0;
dEb1_sum = 0;
dEb2_sum = 0;
for j=1:1:7
for l = 1:1:7
dEa1_sum = dEa1_sum + SO2a1*(w(j,l)*O3_sum - O3w_sum);
dEa2_sum = dEa2_sum + SO2a2*(w(l,j)*O3_sum - O3w_sum);
dEb1_sum = dEb1_sum + SO2b1*(w(j,l)*O3_sum - O3w_sum);
dEb2_sum = dEb2_sum + SO2b2*(w(l,j)*O3_sum - O3w_sum);
end
dE_da(1,j) = -(yd -y(k))*dy_du*2*((O1(1,j)-a(1,j))*O2(1,j))/((b(1,j))*(O3_sum^2))*(dEa1_sum);
dE_da(2,j) = -(yd -y(k))*dy_du*2*((O1(2,j)-a(2,j))*O2(2,j))/((b(2,j)^2*(O3_sum^2)))*(dEa2_sum);
dE_db(1,j) = -(yd -y(k))*dy_du*2*((O1(1,j)-a(1,j))^2*O2(1,j))/((b(1,j)^3*(O3_sum^2)))*(dEb1_sum);
dE_db(2,j) = -(yd -y(k))*dy_du*2*((O1(2,j)-a(2,j))^2*O2(2,j))/((b(2,j)^3*O3_sum^2))*(dEb2_sum);
end
%更新 w a b 三个矩阵（这里将偏差项给省略了）
for i =1:1:7
for j =1:1:7
w(i,j) = w(i,j) - eta*dE_dw(i,j);
end
end
for i = 1:1:2
for j =1:1:7
a(i,j) = a(i,j) - eta*dE_da(i,j);
b(i,j) = b(i,j) - eta*dE_db(i,j);
end
end
end
xlim([0,10]);
figure(1)
plot(e(1,:));
hold on
plot(y(1,:));
legend('e','y'); %e:error y:the output
figure(2)
plot(u(1,:));
时，此时，将四阶龙哥库塔中时间加限制调件，其程序如下：
function dotstate=runge_kutta(~,state)
global u k h;
t=k*h;
dotstate(1:2,1) = zeros(2,1);
Da = 0.072;
phi = 20;
if t <= 5 %时间5分钟之前
B = 8;
delta = 0.3;
else %时间5分钟之后
B = 7.5;
delta = 0.35;
end
dotstate(1) = - state(1) + Da*(1-state(1))*exp(state(2)/(1+state(2)/phi));
dotstate(2) = -(1+delta)*state(2) + B*Da*(1-state(1))*exp(state(2)/(1+state(2)/phi))+delta*u(k);
end
1210
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
