一、理论基础

模糊神经网络（Fuzzy Neural Network，简称FNN）本质上是一种将模糊理论与人工前向神经网络相结合的多层前向神经网络，在处理信息时，该网络能够具有更大的处理范围以及更快的信息处理速度，因此该网络的自学习能力与映射也相对较高。与反向传播神经网络（BPNN）相似，其训练过程主要分为前向传播与反向传播两个部分。
反向传播神经网络（BPNN）原理参考：
反向传播神经网络（BPNN）的实现（Python，附源码及数据集）

1、模糊神经网络结构

模糊神经网络的结构图如下图所示：

2、前向传播过程

数据输入模糊神经网络总共需要经过五层网络结构，最先是输入层，这一层的节点数量为输入数据的特征维度，即当数据的特征维度为n时，输入层的节点数为n。
然后数据从输入层传递至隶属度函数计算层，在这一层需要使用隶属度函数来计算每个节点的隶属度，每一个节点代表一个隶属度函数，这一层的节点个数为输入变量有可能构成的模糊条件的个数。隶属度函数一般选用高斯函数，公式如下：

其中u_ij为节点的输出，x_i为输入数据，c_ij为隶属函数的中心点，b_ij为隶属函数的宽度向量，m为对输入进行模糊分级的个数。
之后数据传递至规则生成层，这一层的节点数为模糊规则数，这一层的输出为每条规则的适用度，一般采用下列公式进行计算：

紧接着在归一化层对规则生成层的节点输出进行归一化处理，这一层的节点与规则生成层相同，使用如下公式进行归一化：

最后数据在输出层进行加权处理获得最终输出，计算的公式为：

其中ω_j为输出层的连接权值，上述公式是基于网络的输出变量维度为1时而设定的权值，当输出变量的维度增加时，对权值也将进行相应的调整。

3、反向传播过程

模糊神经网络同BPNN一样，一般使用反向传播算法如随机梯度下降法（SGD）对隶属函数的中心点、宽度向量以及输出层的连接权值进行优化。
以输出层的权值为例，其更新公式如下：
其中E为误差函数，μ为学习率。
损失函数原理参考：
机器学习基础知识之损失函数
反向传播原理参考：
神经网络之反向传播算法（梯度、误差反向传播算法BP）

4、建模步骤

以使用模糊神经网络进行数据预测为例，可以将模糊神经网络预测模型的建模步骤总结如下：

1. 根据输入数据的相关特征确定模糊神经网络输入层、隶属度函数计算层、规则生成层、归一化层以及输出层的节点数；
2. 选择一种参数初始化方法对模糊神经网络隶属度函数计算层中隶属度函数的中心点、宽度向量以及输出层的连接权值、阈值进行随机初始化；
3. 数据由输入层输入模糊神经网络，依次经过隶属度函数计算层以及规则生成层计算获得输出；
4. 将规则生成层的输出在归一化层进行归一化处理之后经由输出层的加权处理获得网络的最终前向传播输出；
5. 传递至隐含层后经径向基函数对数据进行非线性转换；
6. 选择一种损失函数对网络的前向传播输出以及目标值进行相关计算得到损失值；
7. 以输出层的损失值计算得到输出层连接权值的梯度，选择一种反向传播算法对它们进行调整；
8. 损失值传递至隶属度函数计算层，同样使用相同的反向传播算法对隶属度函数的中心点以及宽度向量进行调整；
9. 获得一个参数得到更新后的模糊神经网络；
10. 在达到最大迭代次数之前，重复步骤3到步骤9，在达到最大迭代次数后，输出隐含层与输出层参数确定的多层前馈神经网络。

参数初始化方法参考：
神经网络基础知识之参数初始化

二、模糊神经网络的实现

以数据预测为例，下面介绍基于Python实现模糊神经网络的过程。
选用某省市的表层土壤重金属元素数据集作为实验数据，该数据集总共96组，随机选择其中的24组作为测试数据集，72组作为训练数据集。选取重金属Ti的含量作为待预测的输出特征，选取重金属Co、Cr、Mg、Pb作为模型的输入特征。

1、训练过程（FNN.py）

#库的导入
import numpy as np
import pandas as pd
#输入数据的导入
df = pd.read_csv("train.csv")
df.columns = ["Co", "Cr", "Mg", "Pb", "Ti"]
Co = df["Co"]
Co = np.array(Co)
Cr = df["Cr"]
Cr = np.array(Cr)
Mg=df["Mg"]
Mg=np.array(Mg)
Pb = df["Pb"]
Pb =np.array(Pb)
Ti = df["Ti"]
Ti = np.array(Ti)
samplein = np.mat([Co,Cr,Mg,Pb])
#数据归一化，将输入数据压缩至0到1之间，便于计算，后续通过反归一化恢复原始值
sampleinminmax = np.array([samplein.min(axis=1).T.tolist()[0],samplein.max(axis=1).T.tolist()[0]]).transpose()#对应最大值最小值
#待预测数据为Ti
sampleout = np.mat([Ti])
sampleoutminmax = np.array([sampleout.min(axis=1).T.tolist()[0],sampleout.max(axis=1).T.tolist()[0]]).transpose()
sampleinnorm = ((np.array(samplein.T)-sampleinminmax.transpose()[0])/(sampleinminmax.transpose()[1]-sampleinminmax.transpose()[0])).transpose()
sampleoutnorm = ((np.array(sampleout.T)-sampleoutminmax.transpose()[0])/(sampleoutminmax.transpose()[1]-sampleoutminmax.transpose()[0])).transpose()
sampleinnorm = sampleinnorm.transpose()
sampleoutnorm = sampleoutnorm.transpose()
S = 2  #模糊分级个数
T = 16  #规则生成层、归一化层节点数
iteration = 500  #迭代训练次数
learningrate = 0.01  #学习率
#c1、b1分别为李隶属度层的中心点与宽度向量、w1为输出层的权值，sampleinnorm.shape[1]代表输入数据的4个类别
c1 = np.random.uniform(low=-1, high=1, size=(S, sampleinnorm.shape[1]))
b1 = np.random.uniform(low=-1, high=1, size=(S, sampleinnorm.shape[1]))
w1 = np.random.uniform(low=0, high=1, size=(1, T))
sampleinnorm = np.mat(sampleinnorm)
#开始训练
for l in range(iteration):
    print("the iteration is :",l+1)
    # 隶属度函数计算层计算输出
    Y = sampleoutnorm.copy()
    u1 = np.zeros((S, sampleinnorm.shape[1]))
    alpha = []
    alpha2 = []
    a1 = []
    for m in range(sampleinnorm.shape[0]):
        for i in range(S):
            for j in range(sampleinnorm.shape[1]):
                u1[i][j] = np.exp((-1) * ((sampleinnorm[m, j] - c1[i][j]) ** 2) / (b1[i][j] ** 2))
        # 规则生成层计算输出
        alpha1 = np.zeros((T, 1))
        for i in range(S):
            a = 1
            if i == 1:
                a = 8
            for p in range(S):
                b = 1
                if p == 1:
                    b = 4
                for q in range(S):
                    c = 1
                    if q == 1:
                        c = 2
                    for k in range(S):
                        alpha1[i * a + p * b + q * c + k] = u1[i][0] * u1[p][1] * u1[q][2] * u1[k][3]
        alpha.append(alpha1)
        # 归一化层计算输出
        alphasum = np.sum(alpha1)
        alpha2.append(alphasum)
        # 输出层计算输出
        a2 = np.dot(w1,alpha1)
        Y[m] = a2/alphasum
    #计算误差
    err = sampleoutnorm - Y
    loss = np.sum(np.abs(err))
    print(" the loss is :",loss)
    #反向传播，分别计算参数w1、c1、b1的误差项
    deltaw1 = np.zeros((1,T))
    deltac1 = np.zeros((S, sampleinnorm.shape[1]))
    deltab1 = np.zeros((S, sampleinnorm.shape[1]))
    for m in range(sampleinnorm.shape[0]):
        changew1 = ((err[m] * alpha[m])/alpha2[m]).transpose()
        changec2 = np.zeros((S, sampleinnorm.shape[1]))
        changeb2 = np.zeros((S, sampleinnorm.shape[1]))
        z=np.zeros((S, sampleinnorm.shape[1]))
        v=np.zeros((S, sampleinnorm.shape[1]))
        for i in range(S):
            for j in range(sampleinnorm.shape[1]):
                z[i][j]=2*(sampleinnorm[m, j] - c1[i][j])/(b1[i][j] ** 2)
                v[i][j]=2*((sampleinnorm[m, j] - c1[i][j])**2)/(b1[i][j]**3)
        h = alpha[m]
        for u in range(T):
            changec1 = (((w1[:,u]-Y[m])*h[u,:])/alpha2[m])*z
            changec2 = changec1+changec2
            changeb1 = (((w1[:,u]-Y[m])*h[u,:])/alpha2[m])*v
            changeb2 = changeb1+changeb2
        deltac1 = deltac1 + changec2
        deltab1 = deltab1 + changeb2
        deltaw1 = deltaw1 + changew1
    #对三个参数进行更新
    c1 = c1 + learningrate * deltac1
    b1 = b1 + learningrate * deltab1
    w1 = w1 + learningrate * deltaw1
print('更新的w1:',w1)
print('更新的b1:',b1)
print('更新的c1:',c1)
#保存训练后的参数
np.save("c1.npy",c1)
np.save("b1.npy",b1)
np.save("w1.npy",w1)

2、测试过程（test.py）

#库的导入
import numpy as np
import pandas as pd
#输入数据的导入，用于测试数据的归一化与返归一化
df = pd.read_csv("train.csv")
df.columns = ["Co", "Cr", "Mg", "Pb", "Ti"]
Co = df["Co"]
Co = np.array(Co)
Cr = df["Cr"]
Cr = np.array(Cr)
Mg=df["Mg"]
Mg=np.array(Mg)
Pb = df["Pb"]
Pb =np.array(Pb)
Ti = df["Ti"]
Ti = np.array(Ti)
samplein = np.mat([Co,Cr,Mg,Pb])
sampleinminmax = np.array([samplein.min(axis=1).T.tolist()[0],samplein.max(axis=1).T.tolist()[0]]).transpose()#对应最大值最小值
sampleout = np.mat([Ti])
sampleoutminmax = np.array([sampleout.min(axis=1).T.tolist()[0],sampleout.max(axis=1).T.tolist()[0]]).transpose()#对应最大值最小值
#导入FNN.py训练好的参数
c1=np.load('c1.npy')
b1=np.load('b1.npy')
w1=np.load('w1.npy')
c1 = np.mat(c1)
b1 = np.mat(b1)
w1 = np.mat(w1)
#测试数据数量
testnum = 24
#测试数据的导入
df = pd.read_csv("test.csv")
df.columns = ["Co", "Cr", "Mg", "Pb", "Ti"]
Co = df["Co"]
Co = np.array(Co)
Cr = df["Cr"]
Cr = np.array(Cr)
Mg=df["Mg"]
Mg=np.array(Mg)
Pb = df["Pb"]
Pb =np.array(Pb)
Ti = df["Ti"]
Ti = np.array(Ti)
input=np.mat([Co,Cr,Mg,Pb])
inputnorm=(np.array(input.T)-sampleinminmax.transpose()[0])/(sampleinminmax.transpose()[1]-sampleinminmax.transpose()[0])
S = 2  #模糊分级个数
T = 16  ##规则生成层、归一化层节点数
#进行预测
out = np.zeros((inputnorm.shape[0],1))
u = np.zeros((S, inputnorm.shape[1]))
# 隶属度函数计算层计算输出
for m in range(inputnorm.shape[0]):
    for i in range(S):
        for j in range(inputnorm.shape[1]):
            u[i][j] = np.exp((-1) * ((inputnorm[m, j] - c1[i,j]) ** 2) / (b1[i,j] ** 2))
    #规则生成层
    alpha = np.zeros((T, 1))
    for i in range(S):
        a = 1
        if i == 1:
            a = 8
        for p in range(S):
            b = 1
            if p == 1:
                b = 4
            for q in range(S):
                c = 1
                if q == 1:
                    c = 2
                for k in range(S):
                    alpha[i * a + p * b + q * c + k] = u[i][0] * u[p][1] * u[q][2] * u[k][3]
    #归一化层
    alphasum = np.sum(alpha)
    alpha1 = np.zeros((T, 1))
    for d in range(T):
        alpha1[d] = alpha[d] / alphasum
    #输出层
    out[m] = np.dot(w1, alpha1)
#对输出结果进行反归一化
out = np.array(out)
diff = sampleoutminmax[:,1]-sampleoutminmax[:,0]
networkout2 = out*diff+sampleoutminmax[0][0]
networkout2 = np.array(networkout2).transpose()
output1=networkout2.flatten()#降成一维数组
output1=output1.tolist()
for i in range(testnum):
    output1[i] = float('%.2f'%output1[i])
print("the prediction is:",output1)
#将输出结果与真实值进行对比，计算误差
output=Ti
print(output.shape)
rmse = (np.sum(np.square(output-output1))/len(output)) ** 0.5
mae = np.sum(np.abs(output-output1))/len(output)
average_loss1=np.sum(np.abs((output-output1)/output))/len(output)
mape="%.2f%%"%(average_loss1*100)
f1 = 0
for m in range(testnum):
    f1 = f1 + np.abs(output[m]-output1[m])/((np.abs(output[m])+np.abs(output1[m]))/2)
f2 = f1 / testnum
smape="%.2f%%"%(f2*100)
print("the MAE is :",mae)
print("the RMSE is :",rmse)
print("the MAPE is :",mape)
print("the SMAPE is :",smape)
#计算预测值与真实值误差与真实值之比的分布
A=0
B=0
C=0
D=0
E=0
for m in range(testnum):
    y1 = np.abs(output[m]-output1[m])/np.abs(output[m])
    if y1 <= 0.1:
        A = A + 1
    elif y1 > 0.1 and y1 <= 0.2:
        B = B + 1
    elif y1 > 0.2 and y1 <= 0.3:
        C = C + 1
    elif y1 > 0.3 and y1 <= 0.4:
        D = D + 1
    else:
        E = E + 1
print("Ratio <= 0.1 :",A)
print("0.1< Ratio <= 0.2 :",B)
print("0.2< Ratio <= 0.3 :",C)
print("0.3< Ratio <= 0.4 :",D)
print("Ratio > 0.4 :",E)

3、测试结果

注：由于每次初始化生成的参数不同，因此对参数设置相同的神经网络进行多次训练和预测，测试结果不会完全一致，此外测试结果的好坏也会受到节点数、学习率、训练次数等参数的影响。

4、参考源码及实验数据集

参考源码及实验数据集