# 附代码-BPPID算法 -python -c -matlab

附代码-BPPID算法 -python -c -matlab

1.算法原理

参考资料：
《先进PID控制matalb仿真.pdf》第二版174页
如下：

其实就是正反馈输出，负反馈根据误差进行权值调整，再次进行计算。

2.相关经验

与传统的 PID 控制方法比较，基于神经网络的 PID 控制器具有如下特点：

(1)由于神经网络可以任意逼近非线性函数，故在选择合适的 BP 网络的隐层数和输出层激活函数的情况下，算法具有很强的逼近性能。 
(2)神经网络具有很强的自学习能力，这体现在网络内部的权值调整上，因此，对 PID 控制参数可以实现在线整定和优化，避免了人工调试 PID 控制参数的繁琐工作，为 PID 控制器的智能化发展提供一条切实可行的途径。 
(3)神经网络具有很强的自适应能力，因此基于神经网络的 PID 控制器对受控对象参数以及结构的变化和输入参考信号的变化具有较强的自适应能力和鲁棒性。

---

当然BPPID也有些许问题，如BP 神经网络 PID 学习速度慢、出现过饱和现象和网络波动震荡等问题，针对这些问题对学习速率、权值修正公式两方面进行优化，可以采用改进学习速率和添加动量修正因子两种方法来解决以上问题。
参考《基于BP神经网络的PID控制算法参数优化的研究_赫煊.pdf》

    初始权值的选取与所选的激活函数有关，例如 tanh 函数，它的输出为**(-1,1)**
    之间的数，则选取的初始权值必须是(-1,1)之间的随机数且均不相同，来减小没有
    必要的重复输出值，可将初始权值靠近变化率比较大的 0 处，来加快网络一开始
    的响应速度，但也不能无限靠近 0，这会使权值变化量Δ𝑤在一开始就过小，从而
    导致训练过慢。

BP 神经网络 PID 算法的缺陷也可通过其他智能控制来弥补，如模糊神经网络 PID 或是粒子群神经网络 PID 等等。

当然下面的代码也会出现目标值修改之后，BP神经网络PID算法无法收敛的问题，需要修改初始矩阵的值来使得BP神经网络PID算法收敛。正确拟合。

3.代码实现

python代码实现：

#PID.py
import numpy as np



class IncrementalPID:
    def __init__(self):

        self.xite_1 = 0.2
        self.alfa = 0.95 #初始化学习率
        self.IN = 4#输入4个值
        self.H = 5#隐含层节点5个
        self.Out = 3 # 输出层3个
        self.wi = np.mat([[-0.6394, -0.2696, -0.3756, -0.7023],
                         [-0.8603, -0.2013, -0.5024, -0.2596],
                         [-1.0000, 0.5543, -1.0000, -0.5437],
                         [-0.3625, -0.0724, 0.6463, -0.2859],
                         [0.1425, 0.0279, -0.5406, -0.7660]]
                         )
        #它是神经网络中输入层到隐含层的权重矩阵。该矩阵是一个5x4的矩阵，包含了20个权重值。这些权重值将用于计算隐含层的输入，是神经网络中的关键参数之一。

        self.wi_1 = self.wi
        self.wi_2 = self.wi
        self.wi_3 = self.wi
        self.wo = np.mat([[0.7576, 0.2616, 0.5820, -0.1416, -0.1325],
                          [-0.1146, 0.2949, 0.8352, 0.2205, 0.4508],
                          [0.7201, 0.4566, 0.7672, 0.4962, 0.3632]]
                         )
        #它是神经网络中隐含层到输出层的权重矩阵。该矩阵是一个3x5的矩阵，包含了15个权重值。

        self.wo_1 = self.wo
        self.wo_2 = self.wo
        self.wo_3 = self.wo

        self.Kp = 0.0
        self.Ki = 0.0
        self.Kd = 0.0

        # self.ts = 40  # 采样周期取值
        self.x = [self.Kp, self.Ki, self.Kd]  # 比例、积分、微分初值

        self.y = 0.0  # 系统输出值
        self.y_1 = 0.0  # 上次系统输出值
        self.y_2 = 0.0  # 上次系统输出值

        self.e = 0.0  # 输出值与输入值的偏差
        self.e_1 = 0.0
        self.e_2 = 0.0
        self.de_1 = 0.0

        self.u = 0.0
        self.u_1 = 0
        self.u_2 = 0.0
        self.u_3 = 0.0
        self.u_4 = 0.0
        self.u_5 = 0.0


        self.Oh = np.mat(np.zeros((self.H, 1)))  # %隐含层的输出  它是一个大小为Hx1的零矩阵。
        self.I = self.Oh  # 隐含层的输入


        self.Oh_sub = [0.0, 0.0, 0.0, 0.0, 0.0]
        self.K_sub = [0.0, 0.0, 0.0]
        self.dK_sub = [0.0, 0.0, 0.0]
        self.delta3_sub = [0.0, 0.0, 0.0]
        self.dO_sub = [0.0, 0.0, 0.0, 0.0, 0.0]
        self.delta2_sub = [0.0, 0.0, 0.0, 0.0, 0.0]

        self.den = -0.8251
        self.num = 0.2099

    # 设置PID控制器参数
    def SetStepSignal(self, SetSignal):

        self.y = -self.den * self.y_1 + self.num * self.u_5
        print(self.y)
        self.e = SetSignal - self.y  #误差值
        print(self.e)
        self.xi = np.mat([SetSignal, self.y, self.e, 5])#PID控制器的输入信号、系统输出值、误差和一个常数。

        self.x[0] = self.e - self.e_1  # 比例输出
        self.x[1] = self.e  # 积分输出
        self.x[2] = self.e - 2 * self.e_1 + self.e_2  # 微分输出
        self.epid = np.mat([[self.x[0]], [self.x[1]], [self.x[2]]])#列 包含self.x[0]、self.x[1]和self.x[2]的列向量

        self.I = np.dot(self.xi, (self.wi.T))#这行代码计算了self.xi和self.wi的转置矩阵的乘积，并将结果赋给了self.I。

        for i1 in range(5):
            self.Oh_sub[i1] = (np.e ** (self.I.tolist()[0][i1]) - np.e ** (-self.I.tolist()[0][i1])) / (np.e ** (self.I.tolist()[0][i1]) + np.e ** (-self.I.tolist()[0][i1]))  # %在激活函数作用下隐含层的输出
        self.Oh = np.mat([[self.Oh_sub[0]], [self.Oh_sub[1]], [self.Oh_sub[2]], [self.Oh_sub[3]], [self.Oh_sub[4]]])

        self.K = np.dot(self.wo, self.Oh)  # 输出层的输入，即隐含层的输出*权值

        for i2 in range(3):
            self.K_sub[i2] = (np.e ** (self.K.tolist()[i2][0])) / (np.e ** (self.K.tolist()[i2][0]) + np.e ** (-self.K.tolist()[i2][0]))  # 输出层的输出，即PID三个参数
        self.K = np.mat([[self.K_sub[0]], [self.K_sub[1]], [self.K_sub[2]]])

        self.Kp = self.K_sub[0]
        self.Ki = self.K_sub[1]
        self.Kd = self.K_sub[2]
        self.Kpid = np.mat([self.Kp, self.Ki, self.Kd]) #行

        self.du = np.dot(self.Kpid, self.epid).tolist()[0][0]
        self.u = self.u_1 + self.du
        # if self.u >= 10:
        #     self.u = 10
        #
        # if self.u <= -10:
        #     self.u = -10
        
        
        # 误差变化量
        self.de = self.e - self.e_1
        if self.de > (self.de_1 * 1.04):
            self.xite = 0.7 * self.xite_1
        elif self.de < self.de_1:
            self.xite = 1.05 * self.xite_1
        else:
            self.xite = self.xite_1

        #  权值在线调整
        self.dyu = np.sin((self.y - self.y_1) / (self.u - self.u_1 + 0.0000001))

        for i3 in range(3):
            self.dK_sub[i3] = 2 / ((np.e ** (self.K_sub[i3]) + np.e ** (-self.K_sub[i3])) * (np.e ** (self.K_sub[i3]) + np.e ** (-self.K_sub[i3])))
        self.dK = np.mat([self.dK_sub[0], self.dK_sub[1], self.dK_sub[2]])


        for i4 in range(3):
            self.delta3_sub[i4] = self.e * self.dyu * self.epid.tolist()[i4][0] * self.dK_sub[i4]
        self.delta3 = np.mat([self.delta3_sub[0], self.delta3_sub[1], self.delta3_sub[2]])

        for l in range(3):
            for i5 in range(5):
                self.d_wo = (1 - self.alfa) * self.xite * self.delta3_sub[l] * self.Oh.tolist()[i5][0] + self.alfa * (self.wo_1 - self.wo_2)

        # self.wo = self.wo_1 + self.d_wo + self.alfa * (self.wo_1 - self.wo_2)
        self.wo = self.wo_1 + self.d_wo

        for i6 in range(5):
            self.dO_sub[i6] = 4 / ((np.e ** (self.I.tolist()[0][i6]) + np.e ** (-self.I.tolist()[0][i6])) * (np.e ** (self.I.tolist()[0][i6]) + np.e ** (-self.I.tolist()[0][i6])))
        self.dO = np.mat([self.dO_sub[0], self.dO_sub[1], self.dO_sub[2], self.dO_sub[3], self.dO_sub[4]])

        self.segma = np.dot(self.delta3, self.wo)

        for i7 in range(5):
             self.delta2_sub[i7] = self.dO_sub[i7] * self.segma.tolist()[0][i7]
        self.delta2 = np.mat([self.delta2_sub[0], self.delta2_sub[1], self.delta2_sub[2], self.delta2_sub[3], self.delta2_sub[4]])

        self.d_wi = (1 - self.alfa) * self.xite * self.delta2.T * self.xi + self.alfa * (self.wi_1 - self.wi_2)
        self.wi = self.wi_1 + self.d_wi

        # 参数更新
        self.u_5 = self.u_4
        self.u_4 = self.u_3
        self.u_3 = self.u_2
        self.u_2 = self.u_1
        self.u_1 = self.u

        self.y_2 = self.y_1
        self.y_1 = self.y

        self.wo_3 = self.wo_2
        self.wo_2 = self.wo_1
        self.wo_1 = self.wo

        self.wi_3 = self.wi_2
        self.wi_2 = self.wi_1
        self.wi_1 = self.wi

        self.e_2 = self.e_1
        self.e_1 = self.e
        self.xite_1 = self.xite

#BPPID.py
import PID
import matplotlib.pyplot as plt

plt.figure(1)  # 创建图表1

def TestPID():
    Process = PID.IncrementalPID()
    
    #画图用
    ProcessXaxis = [0]
    ProcessYaxis = [0]
    X = [0]
    Y = [0]

    for i in range(1, 100):
        print(i)
        #开始调参
        Process.SetStepSignal(10)#设置目标值为15，可以修改
        # Process.SetInertiaTime(100, 50)
        
        #绘图设置
        ProcessXaxis.append(i)
        ProcessYaxis.append(Process.y)
        
        X.append(i)
        Y.append(Process.e)


    plt.figure(1)
    plt.plot(ProcessXaxis, ProcessYaxis, 'r')

    plt.xlim(0, 100)
    plt.axhline(y=10, color='b', linestyle='--')
    plt.ylim(0, 100)
    plt.grid(True)
    plt.title("IncrementalPID")
    plt.show()

if __name__ =='__main__':
    TestPID()

C语言实现代码
C语言实现比较麻烦，需要自己写矩阵库，然后进行矩阵运算

//BPPID.c
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "BPPID.h"

#define H 5 //隐含层参数
#define IN 4 //输入层参数
#define OUT 3 //输出层参数

//创建矩阵函数 行，列，初始化的值 createMatrix(5,4,wi_init)
Matrix createMatrix(int rows, int cols, double initValues[rows][cols]) {
    Matrix mat;
    mat.rows = rows;
    mat.cols = cols;

    mat.data = (double**)malloc(rows * sizeof(double*));
    for (int i = 0; i < rows; i++) {
        mat.data[i] = (double*)malloc(cols * sizeof(double));
        for (int j = 0; j < cols; j++) {
            mat.data[i][j] = initValues[i][j]; // 使用传入的列表值进行初始化
        }
    }

    return mat;
}

//打印矩阵  printMatrix(mat);
void printMatrix(Matrix mat) {
    for (int i = 0; i < mat.rows; i++) {
        for (int j = 0; j < mat.cols; j++) {
            printf("%lf ", mat.data[i][j]);
        }
        printf("\n");
    }
}


//矩阵乘法 matMul(&mat, &trans)
Matrix matMul(Matrix* mat1, Matrix* mat2) {
    Matrix result;
    result.rows = mat1->rows;
    result.cols = mat2->cols;

    result.data = (double**)malloc(result.rows * sizeof(double*));
    for (int i = 0; i < result.rows; i++) {
        result.data[i] = (double*)malloc(result.cols * sizeof(double));
        for (int j = 0; j < result.cols; j++) {
            result.data[i][j] = 0;
            for (int k = 0; k < mat1->cols; k++) {
                result.data[i][j] += mat1->data[i][k] * mat2->data[k][j];
            }
        }
    }

    return result;
}

//矩阵转置 transposeMatrix(&mat)
Matrix transposeMatrix(Matrix* mat) {
    Matrix transposed;
    transposed.rows = mat->cols;
    transposed.cols = mat->rows;

    transposed.data = (double**)malloc(transposed.rows * sizeof(double*));
    for (int i = 0; i < transposed.rows; i++) {
        transposed.data[i] = (double*)malloc(transposed.cols * sizeof(double));
        for (int j = 0; j < transposed.cols; j++) {
            transposed.data[i][j] = mat->data[j][i]; // 转置矩阵
        }
    }

    return transposed;
}
//内存释放
void freeMatrix(Matrix* mat) {
    for (int i = 0; i < mat->rows; i++) {
        free(mat->data[i]);
    }
    free(mat->data);
}


// IncrementalPID 初始化函数
_IncrementalPID Motor_Bppid;


void IncrementalPID_init(void) {
    Motor_Bppid.xite_1 = 0.20;
    Motor_Bppid.alfa = 0.05;

    double wi_init[H][IN] = {
                {-0.2846,0.2193,-0.5097,-1.0668}, 
                {-0.7484,-0.1210,-0.4708,0.0988},
                {-0.7176,0.8297,-1.6000,0.2049},
                {-0.0858,0.1925,-0.6346,0.0347}, 
                {0.4358,0.2369,-0.4564,0.1324}
            };
    Motor_Bppid.wi=createMatrix(5,4,wi_init);//5X4
    Motor_Bppid.wi_1=Motor_Bppid.wi;
    Motor_Bppid.wi_2=Motor_Bppid.wi;
    Motor_Bppid.wi_3=Motor_Bppid.wi;

    double wo_init[OUT][H] = {
                {0.5,0.5, 0.5, 0.5, 0.5},
                {0.5,0.5, 0.5, 0.5, 0.5},
                {0.5,0.5, 0.5, 0.5, 0.5}
            };        
    Motor_Bppid.wo = createMatrix(3,5,wo_init);//3*5
    Motor_Bppid.wo_1=Motor_Bppid.wo;
    Motor_Bppid.wo_2=Motor_Bppid.wo;
    Motor_Bppid.wo_3=Motor_Bppid.wo;

    Motor_Bppid.Kp = 0.0, Motor_Bppid.Ki = 0.0, Motor_Bppid.Kd = 0.0;
    Motor_Bppid.y = 0.0, Motor_Bppid.y_1 = 0.0, Motor_Bppid.y_2 = 0.0;
    Motor_Bppid.e = 0.0, Motor_Bppid.e_1 = 0.0, Motor_Bppid.e_2 = 0.0;
    Motor_Bppid.u = 0.0, Motor_Bppid.u_1 = 0.0, Motor_Bppid.u_2 = 0.0, Motor_Bppid.u_3 = 0.0, Motor_Bppid.u_4 = 0.0, Motor_Bppid.u_5 = 0.0;
    
    double oh_init[1][5] = {0.0, 0.0, 0.0, 0.0, 0.0};
    Motor_Bppid.Oh = createMatrix(1, 5, oh_init);//1*5

    double I_init[1][5] = {0.0, 0.0, 0.0, 0.0, 0.0};
    Motor_Bppid.I = createMatrix(1,5,I_init);

    double oh_sub_init[1][5] = {0.0,0.0,0.0,0.0,0.0};
    Motor_Bppid.Oh_sub = createMatrix(1,5,oh_sub_init);

    double K_sub_init[1][3] = {0.0,0.0,0.0};
    Motor_Bppid.K_sub = createMatrix(1,3,K_sub_init);

    double dK_sub_init[1][3] = {0.0,0.0,0.0};
    Motor_Bppid.dK_sub = createMatrix(1,3,oh_init);
    
    double delta3_sub_init[1][3] = {0.0,0.0,0.0};
    Motor_Bppid.delta3_sub = createMatrix(1,3,delta3_sub_init);

    double dO_sub_init[1][5] = {0.0,0.0,0.0,0.0,0.0};
    Motor_Bppid.dO_sub = createMatrix(1,5,dO_sub_init);

    double delta2_sub_init[1][5] = {0.0,0.0,0.0,0.0,0.0};
    Motor_Bppid.delta2_sub = createMatrix(1,5,delta2_sub_init);

}


// 激活函数
double activationFunction(double x) {
    return (exp(x) - exp(-x)) / (exp(x) + exp(-x));
}
//输出层激活函数
double GactivationFunction(double x) {
    return exp(x) / (exp(x) + exp(-x));
}


// SetStepSignal 函数的C语言版本
double SetStepSignal(double SetSignal, double feedvalue) 
{
    //初始化输入为向量
    // 计算误差值
    Motor_Bppid.e = SetSignal - feedvalue;
    // 示例：初始化矩阵 xi
    double xi_mat_init[1][4] = {SetSignal, feedvalue, Motor_Bppid.e, 5.0};

    Matrix xi_mat= createMatrix(1,4,xi_mat_init);//1*4
    printf("ximat\n");
    printMatrix(xi_mat);

    double epid_mat_init[1][3] = {Motor_Bppid.e - Motor_Bppid.e_1,Motor_Bppid.e,Motor_Bppid.e - 2 * Motor_Bppid.e_1 + Motor_Bppid.e_2};
    Matrix epid_mat= createMatrix(1,3,epid_mat_init);//1*3
    printf("epid_mat\n");
    printMatrix(epid_mat);
    //计算输入层权重
    // 对 wi 进行转置
    Matrix wi_transposed = transposeMatrix(&Motor_Bppid.wi);//4*5
    // 将 xi_mat 与 wi_transposed 相乘
    Matrix result = matMul(&xi_mat, &wi_transposed);// 结果矩阵 1*4 4*5  =  1*5 
    printf(" input weight result\n");
    printMatrix(result);
    //激活
    for (int i = 0; i < H; i++) {
        result.data[0][i] = activationFunction(result.data[0][i]);
    }
    printf(" activate weight result\n");
    printMatrix(result);
    // 将计算结果后的 1*5 矩阵转置为 5*1
    Matrix transposed_result=transposeMatrix(&result); // 转置后的矩阵

    printf(" activate weight transposed result\n");
    printMatrix(transposed_result);
    //计算输出层权重
    Matrix K =matMul(&Motor_Bppid.wo,&transposed_result); // 存储结果的矩阵 // 将 wo 与输出层的转置矩阵相乘，结果存储在 K 矩阵中
    
    printf(" ouput weight result\n");
    printMatrix(K);
    // 调用激活函数对输出层的矩阵进行处理
    // 计算PID参数
    Motor_Bppid.Kp = GactivationFunction(K.data[0][0]);
    Motor_Bppid.Ki = GactivationFunction(K.data[1][0]);
    Motor_Bppid.Kd = GactivationFunction(K.data[2][0]);
    // 将 xi_mat 与 wi 矩阵相乘，结果存储在 result
    printf("Kp: %f\n", Motor_Bppid.Kp);
    // 计算PID增量
    double du =  Motor_Bppid.Kp * epid_mat.data[0][0] + Motor_Bppid.Ki* epid_mat.data[0][1] + Motor_Bppid.Kd * epid_mat.data[0][2];
    printf("du= %f\n", du);
    // 更新输出
    //Motor_Bppid.u = Motor_Bppid.u_1 + du;
    Motor_Bppid.u = feedvalue + du;
    printf("Motor_Bppid.u_1 = %f\n", Motor_Bppid.u_1);
    printf("Motor_Bppid.u = %f\n", Motor_Bppid.u);

    //误差变化量
    // double de = Motor_Bppid.e - Motor_Bppid.e_1;
    // if (de > (de_1 * 1.04)) {
    //     double xite = 0.7 * Motor_Bppid.xite_1;
    // } else if (de < de_1) {
    //     xite = 1.05 * Motor_Bppid.xite_1;
    // } else {
    //     xite = Motor_Bppid.xite_1;
    // }

    // 更新权值调整
    double dyu = sin((feedvalue - Motor_Bppid.y_1) / (du - Motor_Bppid.u_1 + 0.0000001));
    printf("dyu: %f\n", dyu);
    for (int i = 0; i < 3; i++) {
        Motor_Bppid.dK_sub.data[0][i] = 2 / ((exp(Motor_Bppid.K_sub.data[0][i]) + exp(-Motor_Bppid.K_sub.data[0][i])) * (exp(Motor_Bppid.K_sub.data[0][i]) + exp(-Motor_Bppid.K_sub.data[0][i])));
    }
    printf(" Motor_Bppid.dK_sub\n");
    printMatrix(Motor_Bppid.dK_sub);
    // 假设epid已经定义并初始化
    for (int i = 0; i < 3; i++) {
        Motor_Bppid.delta3_sub.data[0][i] = Motor_Bppid.e * dyu * epid_mat.data[0][i] * Motor_Bppid.dK_sub.data[0][i]; // 注意epid的定义和使用
    }

    printf(" Motor_Bppid.delta3_sub\n");
    printMatrix(Motor_Bppid.delta3_sub);
    //输出层
    // 更新wo的权值
    double d_wo_init[3][5] = {
                {0.0, 0.0, 0.0, 0.0, 0.0},
                {0.0, 0.0, 0.0, 0.0, 0.0},
                {0.0, 0.0, 0.0, 0.0, 0.0}
            };   
    Matrix d_wo = createMatrix(3, 5,d_wo_init);


    for (int l = 0; l < 3; l++) {
        for (int i = 0; i < 5; i++) {
            //d_wo.data[l][i] = (1-Motor_Bppid.alfa) * Motor_Bppid.xite_1 * Motor_Bppid.delta3_sub.data[0][l] * result.data[0][i] + Motor_Bppid.alfa * (Motor_Bppid.wo_1.data[l][i] - Motor_Bppid.wo_2.data[l][i]);
            d_wo.data[l][i] =  Motor_Bppid.xite_1 * Motor_Bppid.delta3_sub.data[0][l] * result.data[0][i]+ Motor_Bppid.alfa * (Motor_Bppid.wo_1.data[l][i] - Motor_Bppid.wo_2.data[l][i]);
        }
    }
    printf(" d_wo\n");
    printMatrix(d_wo);
    // 更新wo
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 5; j++) {
            Motor_Bppid.wo.data[i][j] = Motor_Bppid.wo_1.data[i][j] + 0.0001*d_wo.data[i][j]+ Motor_Bppid.alfa * (Motor_Bppid.wo_1.data[i][j] - Motor_Bppid.wo_2.data[i][j]);;
        }
    }
    printf(" change weight result\n");
    printMatrix(Motor_Bppid.wo);
    //输入层
    // 计算dO和delta2
    for (int i = 0; i < 5; i++) {
        Motor_Bppid.dO_sub.data[0][i] = 4 / ((exp(Motor_Bppid.I.data[0][i]) + exp(-Motor_Bppid.I.data[0][i])) * (exp(Motor_Bppid.I.data[0][i]) + exp(-Motor_Bppid.I.data[0][i])));
    }

    printf(" Motor_Bppid.dO_sub\n");
    printMatrix(Motor_Bppid.dO_sub);
    Matrix segma=matMul(&Motor_Bppid.delta3_sub, &Motor_Bppid.wo); // 存储结果的矩阵
    // segma的计算已经完成
    for (int i = 0; i < 5; i++) {
        Motor_Bppid.delta2_sub.data[0][i] = Motor_Bppid.dO_sub.data[0][i] * segma.data[0][i]*0.001; // 注意segma的定义和使用
    }
    printf(" Motor_Bppid.delta2_sub\n");
    printMatrix(Motor_Bppid.delta2_sub);
    // 更新wi的权值
    double d_wi_init[5][4] = {
                {0.0, 0.0, 0.0, 0.0},
                {0.0, 0.0, 0.0, 0.0},
                {0.0, 0.0, 0.0, 0.0},
                {0.0, 0.0, 0.0, 0.0},
                {0.0, 0.0, 0.0, 0.0}
            };   
    //转置delta2
    Matrix delta2_trans = transposeMatrix(&Motor_Bppid.delta2_sub);
    Matrix d_wi = createMatrix(5, 4,d_wi_init);
    for (int i = 0; i < 5; i++) {
        for (int j = 0; j < 4; j++) { // 假设xi是4维的
            d_wi.data[i][j] = Motor_Bppid.xite_1 * delta2_trans.data[i][0] * xi_mat.data[0][j] + Motor_Bppid.alfa * (Motor_Bppid.wi_1.data[i][j] - Motor_Bppid.wi_2.data[i][j]);
        }
    }
    printf(" d_wi\n");
    printMatrix(d_wi);
    // 更新wi
    for (int i = 0; i < 5; i++) {
        for (int j = 0; j < 4; j++) {
            Motor_Bppid.wi.data[i][j] = Motor_Bppid.wi_1.data[i][j] + 0.0001*d_wi.data[i][j];
        }
    }
    printf(" Motor_Bppid.wi\n");
    printMatrix(Motor_Bppid.wi);



    Motor_Bppid.u_5 = Motor_Bppid.u_4;
    Motor_Bppid.u_4 = Motor_Bppid.u_3;
    Motor_Bppid.u_3 = Motor_Bppid.u_2;
    Motor_Bppid.u_2 = Motor_Bppid.u_1;
    Motor_Bppid.u_1 = Motor_Bppid.u;

    Motor_Bppid.y_2 = Motor_Bppid.y_1;
    Motor_Bppid.y_1 = feedvalue;

    Motor_Bppid.wo_3 = Motor_Bppid.wo_2;
    Motor_Bppid.wo_2 = Motor_Bppid.wo_1;
    Motor_Bppid.wo_1 = Motor_Bppid.wo;

    Motor_Bppid.wi_3 = Motor_Bppid.wi_2;
    Motor_Bppid.wi_2 = Motor_Bppid.wi_1;
    Motor_Bppid.wi_1 = Motor_Bppid.wi;

    Motor_Bppid.e_2 = Motor_Bppid.e_1;
    Motor_Bppid.e_1 = Motor_Bppid.e;
    Motor_Bppid.xite_1 = Motor_Bppid.xite_1;
    
    printf("Motor_Bppid.u = %f\n", Motor_Bppid.u);


    return Motor_Bppid.u;

}

int main() {
    //参数初始化
    double value;
    double value_e;
    IncrementalPID_init();
    // double target = 100.0;
    // double value = 2.0;
    // double out_U = SetStepSignal(target,value); 
    // printf("ouput next value %lf: ", out_U);
    for (int i = 0; i < 10; i++) {
        printf("Enter actual value %d: ", i+1);
        scanf("%lf", &value); // 从标准输入读取实际值
        double target = 100.0;
        value_e = value;
        double out_U = SetStepSignal(target,value_e); // 将实际值存储在数组中
        printf("ouput next value %lf: ", out_U);
    }
    

    return 0;
}

matlab代码：

%BP based PID Control 
clear all; 
close all; 
xite=0.20; 
alfa=0.05; 
S=2; %Signal type 
IN=4;H=5;Out=3; %NN Structure 
if S==1 %Step Signal 
wi=[-0.6394 -0.2696 -0.3756 -0.7023; 
 -0.8603 -0.2013 -0.5024 -0.2596; 
 -1.0749 0.5543 -1.6820 -0.5437; 
 -0.3625 -0.0724 -0.6463 -0.2859; 
 0.1425 0.0279 -0.5406 -0.7660]; 
%wi=0.50*rands(H,IN); 
wi_1=wi;wi_2=wi;wi_3=wi; 
wo=[0.7576 0.2616 0.5820 -0.1416 -0.1325; 
 -0.1146 0.2949 0.8352 0.2205 0.4508; 
 0.7201 0.4566 0.7672 0.4962 0.3632]; 
%wo=0.50*rands(Out,H); 
wo_1=wo;wo_2=wo;wo_3=wo; 
end 
if S==2 %Sine Signal 
wi=[-0.2846 0.2193 -0.5097 -1.0668; 
 -0.7484 -0.1210 -0.4708 0.0988; 
 -0.7176 0.8297 -1.6000 0.2049; 
 -0.0858 0.1925 -0.6346 0.0347; 
 0.4358 0.2369 -0.4564 -0.1324]; 
%wi=0.50*rands(H,IN); 
wi_1=wi;wi_2=wi;wi_3=wi; 
wo=[1.0438 0.5478 0.8682 0.1446 0.1537; 
 0.1716 0.5811 1.1214 0.5067 0.7370; 
 1.0063 0.7428 1.0534 0.7824 0.6494]; 
%wo=0.50*rands(Out,H); 
wo_1=wo;wo_2=wo;wo_3=wo; 
end
x=[0,0,0]; 
du_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; 
Oh=zeros(H,1); %Output from NN middle layer 
I=Oh; %Input to NN middle layer 
error_2=0; 
error_1=0; 
ts=0.001; 
for k=1:1:6000 
time(k)=k*ts; 
if S==1 
 rin(k)=1.0; 
elseif S==2 
 rin(k)=sin(1*2*pi*k*ts); 
end 
%Unlinear model 
a(k)=1.2*(1-0.8*exp(-0.1*k)); 
yout(k)=a(k)*y_1/(1+y_1^2)+u_1; 
error(k)=rin(k)-yout(k); 
xi=[rin(k),yout(k),error(k),1]; 
x(1)=error(k)-error_1; 
x(2)=error(k); 
x(3)=error(k)-2*error_1+error_2; 
epid=[x(1);x(2);x(3)]; 
I=xi*wi'; 
for j=1:1:H 
 Oh(j)=(exp(I(j))-exp(-I(j)))/(exp(I(j))+exp(-I(j))); %Middle Layer 
end 
K=wo*Oh; %Output Layer 
for l=1:1:Out 
 K(l)=exp(K(l))/(exp(K(l))+exp(-K(l))); %Getting kp,ki,kd 
end 
kp(k)=K(1);ki(k)=K(2);kd(k)=K(3); 
Kpid=[kp(k),ki(k),kd(k)]; 
du(k)=Kpid*epid; 
u(k)=u_1+du(k); 
dyu(k)=sign((yout(k)-y_1)/(du(k)-du_1+0.0001)); 
%Output layer 
for j=1:1:Out 
 dK(j)=2/(exp(K(j))+exp(-K(j)))^2; 
end 
for l=1:1:Out 
 delta3(l)=error(k)*dyu(k)*epid(l)*dK(l); 
end 
for l=1:1:Out 
 for i=1:1:H 
 d_wo=xite*delta3(l)*Oh(i)+alfa*(wo_1-wo_2); 
 end 
end 
 wo=wo_1+d_wo+alfa*(wo_1-wo_2); 
%Hidden layer 
for i=1:1:H 
 dO(i)=4/(exp(I(i))+exp(-I(i)))^2; 
end 
 segma=delta3*wo; 
for i=1:1:H 
 delta2(i)=dO(i)*segma(i); 
end 
d_wi=xite*delta2'*xi; 
wi=wi_1+d_wi+alfa*(wi_1-wi_2); 
%Parameters Update 
du_1=du(k); 
u_5=u_4;u_4=u_3;u_3=u_2;u_2=u_1;u_1=u(k); 
y_2=y_1;y_1=yout(k); 
 
wo_3=wo_2; 
wo_2=wo_1; 
wo_1=wo; 
wi_3=wi_2; 
wi_2=wi_1; 
wi_1=wi; 
error_2=error_1; 
error_1=error(k); 
end 
figure(1); 
plot(time,rin,'r',time,yout,'b'); 
xlabel('time(s)');ylabel('rin,yout'); 
figure(2); 
plot(time,error,'r'); 
xlabel('time(s)');ylabel('error'); 
figure(3); 
plot(time,u,'r'); 
xlabel('time(s)');ylabel('u'); 
figure(4); 
subplot(311); 
plot(time,kp,'r'); 
xlabel('time(s)');ylabel('kp'); 
subplot(312); 
plot(time,ki,'g'); 
xlabel('time(s)');ylabel('ki'); 
subplot(313); 
plot(time,kd,'b'); 
xlabel('time(s)');ylabel('kd');