精华内容
下载资源
问答
  • 解决了5参数模型系列。 在本文中,我们使用表示理论方法来求解一般的13参数高斯模型,该模型可被视为零维量子场论。 我们用表示理论参数来表示动作中的两个线性和十一个二次项。 根据矩阵变量的适当线性组合,这些...
  • 基于LSTM的股票预测模型_python实现_超详细

    万次阅读 多人点赞 2019-07-05 22:25:13
    文章目录一、背景二、主要技术介绍1、RNN模型2、LSTM模型3、控制门工作原理四、代码实现五、案例分析六、参数设置七、结论 一、背景 近年来,股票预测还处于一个很热门的阶段,因为股票市场的波动十分巨大,随时...

    一、背景

    近年来,股票预测还处于一个很热门的阶段,因为股票市场的波动十分巨大,随时可能因为一些新的政策或者其他原因,进行大幅度的波动,导致自然人股民很难对股票进行投资盈利。因此本文想利用现有的模型与算法,对股票价格进行预测,从而使自然人股民可以自己对股票进行预测。
    理论上,股票价格是可以预测的,但是影响股票价格的因素有很多,而且目前为止,它们对股票的影响还不能清晰定义。这是因为股票预测是高度非线性的,这就要预测模型要能够处理非线性问题,并且,股票具有时间序列的特性,因此适合用循环神经网络,对股票进行预测。
    虽然循环神经网络(RNN),允许信息的持久化,然而,一般的RNN模型对具备长记忆性的时间序列数据刻画能力较弱,在时间序列过长的时候,因为存在梯度消散和梯度爆炸现象RNN训练变得非常困难。Hochreiter 和 Schmidhuber 提出的长短期记忆( Long Short-Term Memory,LSTM)模型在RNN结构的基础上进行了改造,从而解决了RNN模型无法刻画时间序列长记忆性的问题。
    综上所述,深度学习中的LSTM模型能够很好地刻画时间序列的长记忆性。

    二、主要技术介绍

    1、RNN模型

    在传统的RNN(循环神经网络)中,所有的w都是同一个w,经过同一个cell的时候,都会保留输入的记忆,再加上另外一个要预测的输入,所以预测包含了之前所有的记忆加上此次的输入。所有RNN都具有一种重复神经网络模块的链式的形式。在标准的RNN中,这个重复的模块只有一个非常简单的结构,例如一个tanh层。
    当权中大于1时,反向传播误差时,误差将会一直放大,导致梯度爆炸;当权中小于1时,误差将会一直缩小,导致梯度消失,进而导致网络权重更新缓慢,无法体现出RNN的长期记忆的效果,使得RNN太过健忘。RNN模型的结构如图:
    在这里插入图片描述

    2、LSTM模型

    长短期记忆模型(long-short term memory)是一种特殊的RNN模型,是为了解决反向传播过程中存在梯度消失和梯度爆炸现象,通过引入门(gate)机制,解决了RNN模型不具备的长记忆性问题,LSTM模型的结构如图:

    图2
    具体来说,LSTM模型的1个神经元包含了1个细胞状态(cell)和3个门(gate)机制。细胞状态(cell)是LSTM模型的关键所在,类似于存储器,是模型的记忆空间。细胞状态随着时间而变化,记录的信息由门机制决定和更新。门机制是让信息选择式通过的方法,通过sigmoid函数和点乘操作实现。sigmoid取值介于0~1之间,乘即点乘则决定了传送的信息量(每个部分有多少量可以通过),当sigmoid取0时表示舍弃信息,取1时表示完全传输(即完全记住)[2]。
    LSTM 拥有三个门,来保护和控制细胞状态:遗忘门(forget gate)、更新门(update gate)和输出门(output gate)。
    细胞状态类似于传送带。直接在整个链上运行,只有一些少量的线性交互。信息在上面流传保持不变会很容易。
    如图:
    在这里插入图片描述

    3、控制门工作原理

    遗忘门
    在这里插入图片描述
    更新门
    在这里插入图片描述
    在这里插入图片描述

    输出门
    在这里插入图片描述

    四、代码实现

    UI

    demo.py
    import tensorflow as tf
    import numpy as np
    import tkinter as tk
    from tkinter import filedialog
    import time
    import pandas as pd
    
    import stock_predict as pred
    
    
    def creat_windows():
        win = tk.Tk()  # 创建窗口
        sw = win.winfo_screenwidth()
        sh = win.winfo_screenheight()
        ww, wh = 800, 450
        x, y = (sw - ww) / 2, (sh - wh) / 2
        win.geometry("%dx%d+%d+%d" % (ww, wh, x, y - 40))  # 居中放置窗口
    
        win.title('LSTM股票预测')  # 窗口命名
    
        f_open =open('dataset_2.csv')
        canvas = tk.Label(win)
        canvas.pack()
    
        var = tk.StringVar()  # 创建变量文字
        var.set('选择数据集')
        tk.Label(win, textvariable=var, bg='#C1FFC1', font=('宋体', 21), width=20, height=2).pack()
    
        tk.Button(win, text='选择数据集', width=20, height=2, bg='#FF8C00', command=lambda: getdata(var, canvas),
                  font=('圆体', 10)).pack()
    
        canvas = tk.Label(win)
        L1 = tk.Label(win, text="选择你需要的 列(请用空格隔开,从0开始)")
        L1.pack()
        E1 = tk.Entry(win, bd=5)
        E1.pack()
        button1 = tk.Button(win, text="提交", command=lambda: getLable(E1))
        button1.pack()
        canvas.pack()
        win.mainloop()
    
    def getLable(E1):
        string = E1.get()
        print(string)
        gettraindata(string)
    
    def getdata(var, canvas):
        global file_path
        file_path = filedialog.askopenfilename()
        var.set("注,最后一个为label")
        # 读取文件第一行标签
        with open(file_path, 'r', encoding='gb2312') as f:
        # with open(file_path, 'r', encoding='utf-8') as f:
            lines = f.readlines()  # 读取所有行
            data2 = lines[0]
        print()
    
        canvas.configure(text=data2)
        canvas.text = data2
    
    def gettraindata(string):
        f_open = open(file_path)
        df = pd.read_csv(f_open)  # 读入股票数据
        list = string.split()
        print(list)
        x = len(list)
        index=[]
        # data = df.iloc[:, [1,2,3]].values  # 取第3-10列 (2:10从2开始到9)
        for i in range(x):
            q = int(list[i])
            index.append(q)
        global data
        data = df.iloc[:, index].values
        print(data)
        main(data)
    
    def main(data):
        pred.LSTMtest(data)
        var.set("预测的结果是:" + answer)
    
    if __name__ == "__main__":
        creat_windows()
    

    stock_predict.py

    import numpy as np
    import matplotlib.pyplot as plt
    import tensorflow as tf
    import pandas as pd
    import math
    
    def LSTMtest(data):
    
        n1 = len(data[0]) - 1 #因为最后一位为label
        n2 = len(data)
        print(n1, n2)
    
        # 设置常量
        input_size = n1  # 输入神经元个数
        rnn_unit = 10    # LSTM单元(一层神经网络)中的中神经元的个数
        lstm_layers = 7  # LSTM单元个数
        output_size = 1  # 输出神经元个数(预测值)
        lr = 0.0006      # 学习率
    
        train_end_index = math.floor(n2*0.9)  # 向下取整
        print('train_end_index', train_end_index)
        # 前90%数据作为训练集,后10%作为测试集
        # 获取训练集
        # time_step 时间步,batch_size 每一批次训练多少个样例
        def get_train_data(batch_size=60, time_step=20, train_begin=0, train_end=train_end_index):
            batch_index = []
            data_train = data[train_begin:train_end]
            normalized_train_data = (data_train - np.mean(data_train, axis=0)) / np.std(data_train, axis=0)  # 标准化
            train_x, train_y = [], []  # 训练集
            for i in range(len(normalized_train_data) - time_step):
                if i % batch_size == 0:
                    # 开始位置
                    batch_index.append(i)
                    # 一次取time_step行数据
                # x存储输入维度(不包括label) :X(最后一个不取)
                # 标准化(归一化)
                x = normalized_train_data[i:i + time_step, :n1]
                # y存储label
                y = normalized_train_data[i:i + time_step, n1, np.newaxis]
                # np.newaxis分别是在行或列上增加维度
                train_x.append(x.tolist())
                train_y.append(y.tolist())
            # 结束位置
            batch_index.append((len(normalized_train_data) - time_step))
            print('batch_index', batch_index)
            # print('train_x', train_x)
            # print('train_y', train_y)
            return batch_index, train_x, train_y
    
        # 获取测试集
        def get_test_data(time_step=20, test_begin=train_end_index+1):
            data_test = data[test_begin:]
            mean = np.mean(data_test, axis=0)
            std = np.std(data_test, axis=0)  # 矩阵标准差
            # 标准化(归一化)
            normalized_test_data = (data_test - np.mean(data_test, axis=0)) / np.std(data_test, axis=0)
            # " // "表示整数除法。有size个sample
            test_size = (len(normalized_test_data) + time_step - 1) // time_step
            print('test_size$$$$$$$$$$$$$$', test_size)
            test_x, test_y = [], []
            for i in range(test_size - 1):
                x = normalized_test_data[i * time_step:(i + 1) * time_step, :n1]
                y = normalized_test_data[i * time_step:(i + 1) * time_step, n1]
                test_x.append(x.tolist())
                test_y.extend(y)
            test_x.append((normalized_test_data[(i + 1) * time_step:, :n1]).tolist())
            test_y.extend((normalized_test_data[(i + 1) * time_step:, n1]).tolist())
            return mean, std, test_x, test_y
    
        # ——————————————————定义神经网络变量——————————————————
        # 输入层、输出层权重、偏置、dropout参数
        # 随机产生 w,b
        weights = {
            'in': tf.Variable(tf.random_normal([input_size, rnn_unit])),
            'out': tf.Variable(tf.random_normal([rnn_unit, 1]))
        }
        biases = {
            'in': tf.Variable(tf.constant(0.1, shape=[rnn_unit, ])),
            'out': tf.Variable(tf.constant(0.1, shape=[1, ]))
        }
        keep_prob = tf.placeholder(tf.float32, name='keep_prob')  # dropout 防止过拟合
    
        # ——————————————————定义神经网络——————————————————
        def lstmCell():
            # basicLstm单元
            # tf.nn.rnn_cell.BasicLSTMCell(self, num_units, forget_bias=1.0,
            # tate_is_tuple=True, activation=None, reuse=None, name=None) 
            # num_units:int类型,LSTM单元(一层神经网络)中的中神经元的个数,和前馈神经网络中隐含层神经元个数意思相同
            # forget_bias:float类型,偏置增加了忘记门。从CudnnLSTM训练的检查点(checkpoin)恢复时,必须手动设置为0.0。
            # state_is_tuple:如果为True,则接受和返回的状态是c_state和m_state的2-tuple;如果为False,则他们沿着列轴连接。后一种即将被弃用。
            # (LSTM会保留两个state,也就是主线的state(c_state),和分线的state(m_state),会包含在元组(tuple)里边
            # state_is_tuple=True就是判定生成的是否为一个元组)
            #   初始化的 c 和 a 都是zero_state 也就是都为list[]的zero,这是参数state_is_tuple的情况下
            #   初始state,全部为0,慢慢的累加记忆
            # activation:内部状态的激活函数。默认为tanh
            # reuse:布尔类型,描述是否在现有范围中重用变量。如果不为True,并且现有范围已经具有给定变量,则会引发错误。
            # name:String类型,层的名称。具有相同名称的层将共享权重,但为了避免错误,在这种情况下需要reuse=True.
            #
    
            basicLstm = tf.nn.rnn_cell.BasicLSTMCell(rnn_unit, forget_bias=1.0, state_is_tuple=True)
            # dropout 未使用
            drop = tf.nn.rnn_cell.DropoutWrapper(basicLstm, output_keep_prob=keep_prob)
            return basicLstm
    
       
    
        def lstm(X):  # 参数:输入网络批次数目
            batch_size = tf.shape(X)[0]
            time_step = tf.shape(X)[1]
            w_in = weights['in']
            b_in = biases['in']
    
            # 忘记门(输入门)
            # 因为要进行矩阵乘法,所以reshape
            # 需要将tensor转成2维进行计算
            input = tf.reshape(X, [-1, input_size])
            input_rnn = tf.matmul(input, w_in) + b_in
            # 将tensor转成3维,计算后的结果作为忘记门的输入
            input_rnn = tf.reshape(input_rnn, [-1, time_step, rnn_unit])
            print('input_rnn', input_rnn)
            # 更新门
            # 构建多层的lstm
            cell = tf.nn.rnn_cell.MultiRNNCell([lstmCell() for i in range(lstm_layers)])
            init_state = cell.zero_state(batch_size, dtype=tf.float32)
    
            # 输出门
            w_out = weights['out']
            b_out = biases['out']
            # output_rnn是最后一层每个step的输出,final_states是每一层的最后那个step的输出
            output_rnn, final_states = tf.nn.dynamic_rnn(cell, input_rnn, initial_state=init_state, dtype=tf.float32)
            output = tf.reshape(output_rnn, [-1, rnn_unit])
            # 输出值,同时作为下一层输入门的输入
            pred = tf.matmul(output, w_out) + b_out
            return pred, final_states
    
        # ————————————————训练模型————————————————————
    
        def train_lstm(batch_size=60, time_step=20, train_begin=0, train_end=train_end_index):
            # 于是就有了tf.placeholder,
            # 我们每次可以将 一个minibatch传入到x = tf.placeholder(tf.float32,[None,32])上,
            # 下一次传入的x都替换掉上一次传入的x,
            # 这样就对于所有传入的minibatch x就只会产生一个op,
            # 不会产生其他多余的op,进而减少了graph的开销。
    
            X = tf.placeholder(tf.float32, shape=[None, time_step, input_size])
            Y = tf.placeholder(tf.float32, shape=[None, time_step, output_size])
            batch_index, train_x, train_y = get_train_data(batch_size, time_step, train_begin, train_end)
            # 用tf.variable_scope来定义重复利用,LSTM会经常用到
            with tf.variable_scope("sec_lstm"):
                pred, state_ = lstm(X) # pred输出值,state_是每一层的最后那个step的输出
            print('pred,state_', pred, state_)
    
            # 损失函数
            # [-1]——列表从后往前数第一列,即pred为预测值,Y为真实值(Label)
            #tf.reduce_mean 函数用于计算张量tensor沿着指定的数轴(tensor的某一维度)上的的平均值
            loss = tf.reduce_mean(tf.square(tf.reshape(pred, [-1]) - tf.reshape(Y, [-1])))
            # 误差loss反向传播——均方误差损失
            # 本质上是带有动量项的RMSprop,它利用梯度的一阶矩估计和二阶矩估计动态调整每个参数的学习率。
            # Adam的优点主要在于经过偏置校正后,每一次迭代学习率都有个确定范围,使得参数比较平稳.
            train_op = tf.train.AdamOptimizer(lr).minimize(loss)
            saver = tf.train.Saver(tf.global_variables(), max_to_keep=15)
    
            with tf.Session() as sess:
                # 初始化
                sess.run(tf.global_variables_initializer())
                theloss = []
                # 迭代次数
                for i in range(200):
                    for step in range(len(batch_index) - 1):
                        # sess.run(b, feed_dict = replace_dict)
                        state_, loss_ = sess.run([train_op, loss],
                                            feed_dict={X: train_x[batch_index[step]:batch_index[step + 1]],
                                                       Y: train_y[batch_index[step]:batch_index[step + 1]],
                                                       keep_prob: 0.5})
                                            #  使用feed_dict完成矩阵乘法 处理多输入
                                            #  feed_dict的作用是给使用placeholder创建出来的tensor赋值
    
    
                                            #  [batch_index[step]: batch_index[step + 1]]这个区间的X与Y
                                            #  keep_prob的意思是:留下的神经元的概率,如果keep_prob为0的话, 就是让所有的神经元都失活。
                    print("Number of iterations:", i, " loss:", loss_)
                    theloss.append(loss_)
                print("model_save: ", saver.save(sess, 'model_save2\\modle.ckpt'))
                print("The train has finished")
            return theloss
    
        theloss = train_lstm()
    
        # ————————————————预测模型————————————————————
        def prediction(time_step=20):
    
            X = tf.placeholder(tf.float32, shape=[None, time_step, input_size])
            mean, std, test_x, test_y = get_test_data(time_step)
            # 用tf.variable_scope来定义重复利用,LSTM会经常用到
            with tf.variable_scope("sec_lstm", reuse=tf.AUTO_REUSE):
                pred, state_ = lstm(X)
            saver = tf.train.Saver(tf.global_variables())
            with tf.Session() as sess:
                # 参数恢复(读取已存在模型)
                module_file = tf.train.latest_checkpoint('model_save2')
                saver.restore(sess, module_file)
                test_predict = []
                for step in range(len(test_x) - 1):
                    predict = sess.run(pred, feed_dict={X: [test_x[step]], keep_prob: 1})
                    predict = predict.reshape((-1))
                    test_predict.extend(predict)  # 把predict的内容添加到列表
    
                # 相对误差=(测量值-计算值)/计算值×100%
                test_y = np.array(test_y) * std[n1] + mean[n1]
                test_predict = np.array(test_predict) * std[n1] + mean[n1]
                acc = np.average(np.abs(test_predict - test_y[:len(test_predict)]) / test_y[:len(test_predict)])
                print("预测的相对误差:", acc)
    
                print(theloss)
                plt.figure()
                plt.plot(list(range(len(theloss))), theloss, color='b', )
                plt.xlabel('times', fontsize=14)
                plt.ylabel('loss valuet', fontsize=14)
                plt.title('loss-----blue', fontsize=10)
                plt.show()
                # 以折线图表示预测结果
                plt.figure()
                plt.plot(list(range(len(test_predict))), test_predict, color='b', )
                plt.plot(list(range(len(test_y))), test_y, color='r')
                plt.xlabel('time value/day', fontsize=14)
                plt.ylabel('close value/point', fontsize=14)
                plt.title('predict-----blue,real-----red', fontsize=10)
                plt.show()
    
    
    
        prediction()
    

    五、案例分析

    1、数据说明
    本实验分析了两种股票种类,为某单支股票(6109个连续时间点)数据data2上证综合指数前复权日线(6230个连续时间点,1991年到2016年)数据作为data2,分别保存在两个文件中,将两个数据集的最后一列设定为label。前90%数据作为训练集,后10%作为测试集。
    Data1:
    在这里插入图片描述
    Data2:
    在这里插入图片描述
    本次实验所采用的为LSTM模型:
    输入神经元个数 input_size = 选取列数
    输出神经元个数 output_size = 1 (预测值个数)
    学习率 lr = 0.0006
    随机初始化初始化网络权重

    2、数据预处理
    零-均值规范化(z-score标准化):
    标准化值,是讲集合中单个数与集合的均值相减的结果除以集合的标准差得到的标准化的结果,该方法类似与正态分布标准化转换,转换函数公式为:
    在这里插入图片描述
    公式中x为需要被标准化的原始值,μ为均值,σ为标准差,σ不等于0。
    Z分数标准化处理后的值代表原始值和集合均值之间的举例,以标准差为单位计算。该值存在正负值,低于均值均为辅助,反之则为证书,其范围为[-∞,+∞],数据均值为0,方差为1。
    3、损失函数
    损失函数(Loss function)是用来估量网络模型的预测值X与真实值Y的不一致程度,它是一个非负实值函数,通常用 L(Y,f(x))来表示。损失函数越小,模型的鲁棒性就越好。损失函数是经验风险函数的核心部分,也是结构风险函数的重要组成部分。
    本实验采取十分常用的均方误差损失
    平方损失也可以理解为是最小二乘法,一般在回归问题中比较常见,最小二乘法的基本原理是:最优拟合直线是使各点到回归直线的距离和最小的直线,即平方和最。同时在实际应用中,均方误差也经常被用为衡量模型的标准:
    在这里插入图片描述

    4、误差标准
    相对偏差是指某一次测量的绝对偏差占平均值的百分比。
    在这里插入图片描述
    5、可视化UI
    在这里插入图片描述
    在这里插入图片描述

    六、参数设置

    1、输入维度及迭代次数
    在这里插入图片描述
    由表一可见,输入维度越多,网络训练效果越好;迭代次数在100次时,网络已经比较稳定。
    2、忘记偏置
    在这里插入图片描述
    由表二可见,在data1(单支股票)forget_bias适当减小,即忘记部分信息,网络训练效果有些许提高,但在data2(大盘)中,网络训练效果却有所下滑。个人认为,可能是因为对于单支股票来说,近2天的数据相关程度比较小,而在大盘中,因为近2天的数据相关程度比较大,毕竟有多方面因素影响股价。
    3、LSTM单元数
    在这里插入图片描述

    由表三可见,两个数据集中,LSTM单元数增加的情况下时,网络训练效果反而下降,可以看出,其实股票行情在7天内的的相关联程度比在14天内的情况高,但是有可能是因为forget_bias过大。因此,在进行一组实验,调整forget_bias值进行比较。
    在这里插入图片描述
    由表四可以看出,在相同LSTM单元数的情况下,forget_bias较小时,预测效果较好,我们可以看出,在LSTM单元数较大的情况下,forget_bias应选取比较小的,以免记忆太多无效信息。
    在这里插入图片描述
    由表五可以看出,在data1和data2两个数据集中,LSTM单元数较小的情况下,forget_bias比较大时,预测效果较好,记忆更多相关信息。因此LSTM单元数较小的情况下,forget_bias应选取比较大的,记忆更多相关信息。

    4、可视化结果
    选取数据集data1,迭代次数为200次

    (1)、忘记偏置=1.0 , LSTM单元数 = 2
    在这里插入图片描述

    (2)、忘记偏置=0.7 , LSTM单元数 = 2(表现最好)
    在这里插入图片描述
    (3)、忘记偏置=1.0 , LSTM单元数 = 7
    在这里插入图片描述
    (4)、忘记偏置=1.0 , LSTM单元数 = 14
    在这里插入图片描述

    (5)、忘记偏置=0.7, LSTM单元数 = 7
    在这里插入图片描述
    (6)、忘记偏置=0.4 , LSTM单元数 = 7
    在这里插入图片描述
    (7)、忘记偏置=0.4 , LSTM单元数 = 14
    在这里插入图片描述

    七、结论

    针对以上实验,可以得知,在LSTM模型下的对股票收盘价预测值较为准确和稳定。对LSTM模型进行参数调整,发现迭代次数在100次后,网络模型趋于稳定,说明其是一个较轻量级的网络;在LSTM单元数较大的情况下,forget_bias应选取比较小的,以免记忆太多无效信息;LSTM单元数较小的情况下,forget_bias应选取比较大的,记忆更多相关信息。当然,这和本身的数据集有关。就股票数据集来说,本实验中表现的最优秀的是,忘记偏置为0.7,LSTM神经单元数取2时,网络预测效果最好,说明,在2天内股票序列是比较有价值的,与最后预测值有一定程度的联系。

    完整程序下载

    注意,代码下载后仍需自行调试~
    积分值为5(如果有变为csdn自行修改)——完整代码其实和上面贴出的差别不大,酌情下载~

    https://download.csdn.net/download/zxm_jimin/12126063

    参考文献
    [1] 陈卫华. 基于深度学习的上证综指波动率预测效果比较研究[D].统计与信息论坛,2018.
    [2] Hochreiter & Schmidhuber. Long short-term memory[ J]. Neural Computation, 1997, 9( 8).

    参考博客
    https://blog.csdn.net/jiaoyangwm/article/details/79725445
    https://blog.csdn.net/mylove0414/article/details/56969181

    感谢各位大佬看到最后~
    本文为原创。转载请注明出处。
    注:原理部分,参考了一些文章,博客,如有侵权请联系,我附上原出处。

    展开全文
  • transfer dense121迁移学习模型模板参数不变只训练classifier这一层的参数 ,迁移学习的一个案例模板
  • 将连续的PID控制转换为数字式时,微分环节被用差分代替,积分环节被累加和代替,比例环节则保持不变。差分的实现非常简单,只需要用 e ( k + 1 ) − e ( k ) e(k+1)-e(k) e ( k + 1 ) − e ( k ) 即 e ( k ) − e 1 ...

    0.符号说明

    1. y(k)——系统响应输出的离散值
    2. u(k)——数字PID控制输出的离散值
    3. r(k)——期望输出的离散值(事先已知),在本例中为常数(即阶跃输入)
    4. e(k)——e(k)=r(k)-y(k),为期望值-实际值,是单位负反馈的误差比较信号
      图片来源于百度百科
      注:图片来源于百度百科

    1.如何根据连续系统建立差分方程

    1.1.获取连续系统的传递函数

    线性定常系统的控制中,PID是个非常常见的控制方式,如果可以通过Matlab仿真出PID的控制效果图,那么对系统设计时的实时调试将会容易得多。在这里我们将会以一个利用系统辨识参数的PID设计为为例展示Matlab仿真PID的过程。
    首先需要对一个未知的系统的参数进行辨识,以延迟环节可以忽略不计的电机调速系统为例。将时间戳导入xdata向量,对应的时刻转速导入ydata向量,进行系统辨识

    链接:Matlab的系统辨识

    我们就以上文链接中辨识的系统传递函数为例:
    G ( s ) = 0.998 0.021 s + 1 G(s)=\frac{0.998}{0.021s+1} G(s)=0.021s+10.998因此通过tf函数建立系统结构体如下:

    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
    

    1.2.获取离散系统的传递函数

    由于是数字PID仿真,我们需要选取一个采样时间,本案例选用的是0.005s(注意,采样周期应该小于系统纯滞后时间的0.1倍)。在对其进行数字PID控制前,我们需要将这个系统离散化:

    ts=0.005;  %采样时间=0.005s
    dsys=c2d(sys,ts,'z');      %离散化
    

    dsys即我们根据采样周期离散化的Z变换系统。首先我们需要提取这个Z变化d那系统的参数方便后面的计算:

    [num,den]=tfdata(dsys,'v');%'v'代表强制以向量的格式(默认为元胞数组)输出num和den
    

    1.3.转换为差分方程

    求解出的Z变换表达式为 d s y s = n u m ( 1 ) ⋅ z + n u m ( 2 ) d e n ( 1 ) ⋅ z + d e n ( 2 ) = 0.2114 z − 0.7881 dsys=\frac{num(1)\cdot z +num(2)}{den(1)\cdot z+den(2)}=\frac{0.2114}{z-0.7881} dsys=den(1)z+den(2)num(1)z+num(2)=z0.78810.2114
    在PID仿真的过程中我们需要求解出时域表达式 ,因此需要借助差分方程解决,对于以下的Z变换:

    \begin{equation}
    Y(z)=dsys\cdot U(z)=\frac{num(2)}{den(1)\cdot z+den(2)}\cdot U(z)
    \label{eq:Sample1}
    \end{equation}

    \begin{equation}
    zY(z)+den(2)Y(z)=num(1)zU(z)+num(2)U(z)
    \label{eq:Sample2}
    \end{equation}
    对上式进行反Z变换,可以得到以下的差分方程:

    \begin{equation}
    y(k+1)+den(2)y(k)=num(1)u(k+1)+num(2)u(k)
    \label{eq:Sample3}
    \end{equation}

    \begin{equation}
    y(k+1)=-den(2)y(k)+num(1)u(k+1)+num(2)u(k)
    \label{eq:Sample4}
    \end{equation}
    位置型PID仿真时实际上可以不需要保存前一个数据(u(k)和y(k)),增量型PID必须要保存前一个数据。这里我们使用了位置型PID,但仍然利用 u 1 u_1 u1 y 1 y_1 y1保存了上一个数据,仅仅是为了演示这一过程。\begin{equation}
    y(k+1)=-den(2)y(k)+num(1)u(k+1)+num(2)u(k)
    \end{equation}
    可以转换为下面的式子:
    \begin{equation}
    y(k)=-den(2)y_1+num(1)u(k)+num(2)u_1
    \label{eq:Sample5}
    \end{equation}
    我们的差分方程就这样建立完毕。注意,此差分方程仅仅是描述系统模型的运算规律的,和我们的控制无关。因此是y(k)和u(k)的映射关系。我们下面的控制则是利用负反馈信号e(k)导出u(k)的输出,求解的是控制器u(k)的序列值。

    2.基本PID控制原理

    以位置型PID控制为例。将连续的PID控制转换为数字式时,微分环节被用差分代替,积分环节被累加和代替,比例环节则保持不变。差分的实现非常简单,只需要用 e ( k + 1 ) − e ( k ) e(k+1)-e(k) e(k+1)e(k) e ( k ) − e 1 e(k)-e_1 e(k)e1等效即可。积分的实现在每一次运算的后面都累加原来的误差,即Ee=Ee+e_1;即可。PID的控制器输出 u ( k ) = K p ⋅ e ( k ) + K d ⋅ ( e ( k ) − e 1 ) + K i ⋅ E e u(k)=Kp\cdot e(k)+Kd\cdot (e(k)-e_1)+Ki\cdot Ee u(k)=Kpe(k)+Kd(e(k)e1)+KiEe
    PID控制器构造完毕,我们需要通过r(k)和y(k)得到e(k),再通过e(k)得出u(k),进而再求解出y(k),再结合r(k)求解出e(k),…以此循环,求解出离散的响应点。
    详细的代码如下:

    ts=0.005;  %采样时间=0.005s
    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
    dsys=c2d(sys,ts,'z');      %离散化
    [num,den]=tfdata(dsys,'v');   %
    e_1=0;      %前一时刻的偏差      
    Ee=0;       %累积偏差
    u_1=0.0;    %前一时刻的控制量
    y_1=0;       %前一时刻的输出
    %PID参数
    kp=0.22;    
    ki=0.13;
    kd=0;
    u=zeros(1,1000);%预先分配内存
    time=zeros(1,1000);%时刻点(设定1000个)
    for k=1:1:1000
        time(k)=k*ts;   %时间参数
        r(k)=1500;      %期望值
        y(k)=-1*den(2)*y_1+num(2)*u_1+num(1)*u(k);%系统响应输出序列
        e(k)=r(k)-y(k);   %误差信号
        u(k)=kp*e(k)+ki*Ee+kd*(e(k)-e_1); %系统PID控制器输出序列
        Ee=Ee+e(k);    %误差的累加和
        u_1=u(k);    	%前一个的控制器输出值
        y_1=y(k);    	%前一个的系统响应输出值
        e_1=e(k);		%前一个误差信号的值
    end
    %(仅绘制过渡过程的曲线,x坐标限制为[0,1])
    p1=plot(time,r,'-.');xlim([0,1]);hold on;%指令信号的曲线(即期望输入)
    p2=plot(time,y,'--');xlim([0,1]);%不含积分分离的PID曲线
    hold on;
    

    输出的PID控制曲线如下:
    PID控制

    3.比较PID输出,分析参数产生的影响

    一个基本的PID就完成了。下面如果我们想要知道修改PID的三个参数kp,ki,kd会带来什么效果,只需要在程序中修改即可。为了方便起见,我们建立一个PID的数组,kp,ki,kd每次都取数组的一个值,然后设定一个大循环开始循环仿真。再利用subplot输出子图的方式将所有的PID效果都输出到一个图进行对比。该代码根据上述代码修改已经很容易,PID比较图的代码如下:

    close all
    PID=[0.22,0.13,0;
        0.4,0.13,0;
        0.4,0.25,0;
        0.8,0.23,0.4;
        0.8,0.2,1;
        0.7,0.2,0.9];%初始化PID参数
    for pid=1:1:6
    ts=0.005;  %采样时间=0.005s
    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
    dsys=c2d(sys,ts,'z');      %离散化
    [num,den]=tfdata(dsys,'v');   %
    e_1=0;      %前一时刻的偏差      
    Ee=0;       %累积偏差
    u_1=0.0;    %前一时刻的控制量
    y_1=0;       %前一时刻的输出
    %PID参数
    kp=PID(pid,1);    
    ki=PID(pid,2);
    kd=PID(pid,3);
    u=zeros(1,1000);
    time=zeros(1,1000);
    for k=1:1:1000
        time(k)=k*ts;   %时间参数
        r(k)=1500;      %给定量
        y(k)=-1*den(2)*y_1+num(2)*u_1+num(1)*u(k);
        e(k)=r(k)-y(k);   %偏差
        u(k)=kp*e(k)+ki*Ee+kd*(e(k)-e_1);   
        Ee=Ee+e(k);    
        u_1=u(k);    
        y_1=y(k);    
        e_1=e(k);
    end
    subplot(2,3,pid);
    p1=plot(time,r,'-.');xlim([0,1]);hold on;
    p2=plot(time,y,'--');xlim([0,1]);
    title(['Kp=',num2str(kp),' Ki=',num2str(ki),' Kd= ',num2str(kd)]);
    hold on;
    end
    

    输出的子图矩阵如下:
    PID子图矩阵
    可以发现,修改Kp会造成上升时间的缩短,但是有可能也会带来较大的超调。积分的增加是一个严重的滞后环节,会减小相位裕度,也会带来超调(超调量并不是绝对的,相对于较小的Kp可能会产生较大的超调,而Kp较大时超调会减小(例如第一行的1图和2图的对比))。然而积分的引入也是必要的,否则将会很长时间无法削弱误差e(k)(例如第二行第二个图)。微分的引入相当于一个超前校正,会减少超调,但是过渡的微分很可能会造成尾部振荡,系统逐渐变得不稳定。因此微分和积分之间需要一个平衡,当满足这个平衡的时候,系统几乎没有振荡,同时响应速度也较快。(第一行的图3是积分过多,产生超调,第二行的图1和图3就比较理想)
    综合上述,PID的调节经验可以归结为以下几点:

    • Kp较小时,系统对微分和积分环节的引入较为敏感,积分会引起超调,微分可能会引起振荡,而振荡剧烈的时候超铁也会增加。
    • Kp增大时,积分环节由于滞后产生的超调逐渐减小,此时如果想要继续减少超调可以适当引入微分环节。继续增大Kp系统可能会不太稳定,因此在增加Kp的同时引入Kd减小超调,可以保证在Kp不是很大的情况下也能取得较好的稳态特性和动态性能。
    • Kp较小时,积分环节不宜过大,Kp较大时积分环节也不宜过小(否则调节时间会非常地长),在下面这个例子中我们还会介绍到,当使用分段PID,在恰当的条件下分离积分,可以取得更好的控制效果。原因在于在稳态误差即将满足要求时,消除了系统的滞后。因此系统超调会明显减少。本例中采样的抗积分饱和的方法是遇限削弱积分法。

    4.改进PID算法(遇限削弱积分法)

    遇限削弱积分法的原理是
    u ( k ) > u m a x u(k)>u_{max} u(k)>umax时,若e(k)>0即输出值还未到达指定值,则认为积分会带来滞后,不再积分。
    u ( k ) < 0 u(k)<0 u(k)<0时,若e(k)<0即输出值超过了指定值,则认为积分会带来滞后,不再积分。
    在本案例中认为 u m a x = r ( k ) u_{max}=r(k) umax=r(k)
    改进PID算法如下(需要些两个循环,当然也可以用一个循环,将其中的PID设为一个子过程调用):

    close all
    ts=0.005;  %采样时间=0.005s
    sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
    dsys=c2d(sys,ts,'z');      %离散化
    [num,den]=tfdata(dsys,'v');   %
    e_1=0;      %前一时刻的偏差      
    Ee=0;       %累积偏差
    u_1=0.0;    %前一时刻的控制量
    y_1=0;       %前一时刻的输出
    %PID参数
    kp=0.22;    
    ki=0.13;
    kd=0;
    u=zeros(1,1000);
    time=zeros(1,1000);
    for k=1:1:1000
        time(k)=k*ts;   %时间参数
        r(k)=1500;      %给定量
        y(k)=-1*den(2)*y_1+num(2)*u_1+num(1)*u(k);
        e(k)=r(k)-y(k);   %偏差
        u(k)=kp*e(k)+ki*Ee+kd*(e(k)-e_1);   
        Ee=Ee+e(k);    
        u_1=u(k);    
        y_1=y(k);    
        e_1=e(k);
    end
    p1=plot(time,r,'-.');xlim([0,1]);hold on;
    p2=plot(time,y,'--');xlim([0,1]);
    hold on;
    a=1;%控制积分分离的二值数
    e_1=0;Ee=0;u_1=0.0;y_1=0;%重新初始化       
    for k=1:1:1000
        time(k)=k*ts;   %时间参数
        r(k)=1500;      %给定量
        y(k)=-1*den(2)*y_1+num(2)*u_1;
        e(k)=r(k)-y(k);   %偏差
        u(k)=kp*e(k)+ki*Ee+kd*(e(k)-e_1);   
         if ((u(k)>r(k)) && (e(k)>0))||((u(k)<0) && (e(k)<0))
             a=0;
         else 
             a=1;
         end     
        Ee=Ee+a*e(k);    
        u_1=u(k);    
        y_1=y(k);    
        e_1=e(k);
    end
    p3=plot(time,y,'-');xlim([0,1]);
    title('含积分分离与不含积分分离的对比');
    legend([p1,p2,p3],'指令信号','不含积分分离','含积分分离');
    

    输出的曲线对比图如下:
    积分分离之后的改进PID
    可以发现,系统的超调量明显减少了,调节时间也减少了一点。原因在于我们采用了分段PID的手段,既消除了稳态误差还削弱了积分环节带来的滞后影响。

    5.simulink仿真

    需要的模块名称(不区分大小写)如下:

    • gain(参数分别为0.22和0.13/0.005)
    • sum(参数分别为"|±"和"|++")
    • integrator
    • scope
      注意:本文使用的是离散PID仿真,而simulink使用的是连续系统仿真,转换PID参数时P参数不变,I参数应该除以仿真间隔Ts=0.005,D参数应该乘Ts。

    以表中第一组PI参数为例:
    在这里插入图片描述
    得到的示波器曲线如下:
    在这里插入图片描述

    希望本文对您有帮助,谢谢阅读。

    展开全文
  • 我们在轻子质量和混合角的模块化不变模型中探索带电轻子扇区的替代描述。 除模量外,我们模型的对称破坏扇区还包括普通黄酮。 中微子质量项仅取决于模量,并经过调整以最小化自由参数的数量。 带电的轻子Yukawa...
  • 背景 \qquad 调整PID参数时,我们常常需要按照各种方法查各种表,还要根据建立的数学模型和之前的结果做人工调整。这个过程常常是漫长的,而且最终得到的很可能也不是最优的参数。若系统的参数有改动,又要人工再调...


    阅读前必看:

    1. 本代码基于MATLAB2017a版本,如果版本不同可能会报错
    2. 请从set_para.m文件开始运行,其他M文件(+下载的资源包里面的slx文件)放在一个文件夹
    3. 每次迭代都会打印出来,如果运行时间过长怀疑程序死机可以观察迭代次数是否有变化
    4. 最后会输出3幅图(figure)
    5. slx文件下载链接:https://download.csdn.net/download/weixin_44044411/11609837

    混合粒子群算法链接:https://blog.csdn.net/weixin_44044411/article/details/103638611?spm=1001.2014.3001.5502
    MATALB的自带的粒子群优化器来了:particleswarm

    0.背景

    \qquad 调整PID参数时,我们常常需要按照各种方法查各种表,还要根据建立的数学模型和之前的结果做人工调整。这个过程常常是漫长的,而且最终得到的很可能也不是最优的参数。若系统的参数有改动,又要人工再调一次。而3个参数6个方向的调整,如果输入的变量很多,那么调节的可能性组合就更多了,这个工作常常不是非专业人士所能胜任的。
    \qquad 现在我给大家介绍一种不需要知道内部原理也可以调整PID参数的方法,而且调整的思路是所需即所得——你需要的参数是什么样子的,这个算法就会往这个方向调整PID参数。如果这个系统是多输入多输出的,我需要让A输出的超调小一些,而B输出的调节时间可以尽量大一些,只需要调整我们的评价函数就可以让算法自动调参以满足客户需求了。那这样的算法使用起来是不是很容易呢?

    1.粒子群算法

    1.1.算法简介

    \qquad 粒子群算法(PSO)于1995年提出,和遗传算法一样,也是一种群体迭代算法,不同的是,粒子群算法需要整定的参数更少,不存在交叉和变异过程,所以收敛速度更快。
    \qquad 由于是群体迭代算法,因此该算法一般做成并行的。如果做成串行的,也是可以运行的,就是会随着种群规模的增大越来越慢。对于Python运用熟练的朋友可以利用Python做成并行的,对于不熟悉并行算法的朋友,直接使用MATLAB就好了,因为MATLAB的矩阵运算已经集成了并行算法。
    \qquad 为了让大家理解,我们做一个很形象的比喻。将我们的粒子群比作一群可以互相通讯的鸟,每只鸟不知道食物的具体位置,但是可以根据嗅觉知道离食物的距离和自己在地图中的位置。种群中的任意2只鸟均可以知道对方离食物的距离和对方的具体位置。假设每只鸟的飞行都是有惯性的,而每只鸟都会根据自己距离食物的距离和位置及其他鸟距离食物的距离和位置选择下一次的飞行方向。最终,所有鸟都会在信息交换的基础上,汇聚在食物地点。
    \qquad 这只是对粒子群算法的一个形象的比喻,不必过多计较该模型的现实性,闲言少叙,让我们进入 算 法 步 骤 \color{red}算法步骤

    1.2.算法步骤

    【step1】确定参数维度 N \color{blue}{N} N、参数范围(即每只鸟的初始位置),确定惯性系数 c 1 , c 2 , w c_1,c_2,w c1,c2,w,确定种群规模m,迭代次数n。
    每个粒子的初始位置是随机的,设输入向量 x = ( x 1 , x 2 , . . . , x N ) T x=(x_1,x_2,...,x_N)^T x=(x1,x2,...,xN)T必须满足参数范围限制为:
    { x m i n ( 1 ) < x 1 < x m a x ( 1 ) x m i n ( 2 ) < x 2 < x m a x ( 2 ) . . . x m i n ( N ) < x 1 < x m a x ( N ) \begin{cases}x_{min}^{(1)}<x_1<x_{max}^{(1)} \\x_{min}^{(2)}<x_2<x_{max}^{(2)} \\... \\x_{min}^{(N)}<x_1<x_{max}^{(N)} \end{cases} xmin(1)<x1<xmax(1)xmin(2)<x2<xmax(2)...xmin(N)<x1<xmax(N)
    X m i n = ( x m i n ( 1 ) , x m i n ( 2 ) , . . . x m i n ( N ) ) , X m a x = ( x m a x ( 1 ) , x m a x ( 2 ) , . . . x m a x ( N ) ) \color{blue}X_{min}=(x_{min}^{(1)},x_{min}^{(2)},...x_{min}^{(N)}),X_{max}=(x_{max}^{(1)},x_{max}^{(2)},...x_{max}^{(N)}) Xmin=(xmin(1),xmin(2),...xmin(N)),Xmax=(xmax(1),xmax(2),...xmax(N))
    则输入向量 x x x应满足 X m i n < x < X m a x X_{min}<x<X_{max} Xmin<x<Xmax
    每个粒子的初速度定为0,即 v 0 = 0 v_0=0 v0=0,第 j j j个粒子( 1 ≤ j ≤ m 1≤j≤m 1jm)的下一次迭代的速度 v ( j ) v^{(j)} v(j)由三部分组成:
    v ( j ) = w ⋅ v 0 + c 1 ⋅ r a n d ⋅ ( P ( j ) − X ( j ) ) + c 2 ⋅ r a n d ⋅ ( P G − X ( j ) ) v^{(j)}=w\cdot v_0+c_1\cdot rand\cdot (P^{(j)}-X^{(j)})+c_2\cdot rand\cdot(P_G-X^{(j)}) v(j)=wv0+c1rand(P(j)X(j))+c2rand(PGX(j))
    注:rand是(0,1)的随机数, v 0 v_0 v0代表上一次粒子的速度。
    第一部分为自身惯性因子,因为下一次的迭代次数保留了上一次的速度信息;
    第二个部分为自身最优因子, P ( j ) \color{blue}{P^{(j)}} P(j)为第 j j j个因子在之前所有迭代次数中自适应度最高的位置,可以理解为历史上自身离食物最近的位置。
    第三部分为社会因子, P G \color{blue}{P_G} PG为种群在之前所有迭代次数中自适应度最高的位置,可以理解为历史上所有粒子离食物最近的位置中的最近的一个。
    一般情况下,取 w = 1 , c 1 = c 2 = 2 w=1,c_1=c_2=2 w=1,c1=c2=2,当种群规模较大时,可以让 w w w的值随迭代次数减小以增加收敛速度。
    【step2】按照step1的速度公式计算下一次的速度,并计算下一次的粒子位置。对于第 j j j个粒子,第 k + 1 k+1 k+1次迭代(第 k + 1 k+1 k+1代)的位置 X k + 1 ( j ) \color{blue}{X_{k+1}^{(j)}} Xk+1(j)与第 k k k次迭代的位置 X K ( j ) \color{blue}{X_K^{(j)}} XK(j)与速度 v k ( k + 1 ) \color{blue}{v_k^{(k+1)}} vk(k+1)关系为:
    X k + 1 ( j ) = X k ( j ) + v k ( j + 1 ) ⋅ d t X^{(j)}_{k+1}=X^{(j)}_{k}+v_k^{(j+1)}\cdot dt Xk+1(j)=Xk(j)+vk(j+1)dt
    其中 d t \color{blue}{dt} dt是仿真间隔,一般情况下可以取1,如果希望仿真得慢一点,搜索平滑一点,可以适当减小 d t \color{blue}{dt} dt
    【step3】计算每个粒子的自适应度 F k + 1 ( j ) \color{blue}{F^{(j)}_{k+1}} Fk+1(j),为了计算出step1公式中的 P ( j ) 和 P G \color{blue}{P^{(j)}和P_G} P(j)PG。为了方便起见,我们记前k次计算得到了的 P G P_G PG P G ( k ) P_G^{(k)} PG(k),则第k+1次迭代中我们将适应度最高的粒子位置记为 P G ( k + 1 ) P_G^{(k+1)} PG(k+1),则最终的 P G P_G PG为:
    P G = { P G ( k ) , F ( P G ( k ) ) > F ( P G ( k + 1 ) ) P G ( k + 1 ) , F ( P G ( k ) ) < F ( P G ( k + 1 ) ) P_G=\begin{cases}P_G^{(k)},\qquad F(P_G^{(k)})>F(P_G^{(k+1)}) \\[2ex]P_G^{(k+1)},\quad F(P_G^{(k)})<F(P_G^{(k+1)}) \end{cases} PG=PG(k)F(PG(k))>F(PG(k+1))PG(k+1)F(PG(k))<F(PG(k+1))
    同样,我们记前k次计算得到的第 j j j个粒子的位置为 P k ( j ) P^{(j)}_{k} Pk(j),第k+1次计算得到的第 j j j个粒子的位置为 P k + 1 ( j ) P^{(j)}_{k+1} Pk+1(j),则最终的第 j j j个粒子的历史最优解 P ( j ) P^{(j)} P(j)为:
    P ( j ) = { P k ( j ) , F ( P k ( j ) ) > F ( P k + 1 ( j ) ) P k + 1 ( j ) , F ( P k ( j ) ) < F ( P k + 1 ( j ) ) P^{(j)}=\begin{cases}P_{k}^{(j)},\quad F(P_{k}^{(j)})>F(P_{k+1}^{(j)}) \\[2ex]P_{k+1}^{(j)},\quad F(P_{k}^{(j)})<F(P_{k+1}^{(j)}) \end{cases} P(j)=Pk(j)F(Pk(j))>F(Pk+1(j))Pk+1(j)F(Pk(j))<F(Pk+1(j))
    【step4】更新每个粒子的信息。
    由于我们在step2的位置迭代中已经更新过粒子的位置信息,在step1中的速度迭代中已经更新过粒子的速度信息,而在step3中又更新了每个粒子的历史最优位置 P ( j ) P^{(j)} P(j)及种群最优位置 P G P_G PG,因此实际上如果仅需要知道最优解的情况下我们不需要做这一步。
    但如果需要作出参数随迭代次数的改变的图的话,每次迭代产生最优解 P G ( k ) P_G^{(k)} PG(k)及最优适应度 F ( P G ( k ) ) F(P_G^{(k)}) F(PG(k))需要用数组保存下来。

    1.3.算法举例

    \qquad 根据以上的算法步骤,本人编写了一段粒子群算法的MATLAB代码,这个代码很简单,是无约束问题下的全局粒子群算法。全局粒子群算法收敛速度快,占用内存小,但是容易陷入局部最优解,每次迭代的最优解可能存在细微差异,当然如果这个无约束问题本身是单解的话,这个差异会随着种群规模的扩大和迭代次数增加而逐渐消除。

    function [Pg,fmin]=PSO(f,dimension,n,m,xmax,xmin,vmax,vmin)
    %全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
    %位置限幅为[xmin,xmax],速度限幅为[vmin,vmax]
        Savef=zeros(n+1,1);
        SaveData=zeros(m,dimension,10);%记录11代种群的位置
        w=1;c1=2;c2=2;%速度惯性系数
        dt=0.3;%位移仿真间隔
        v=zeros(m,dimension);%初始速度为0
        X=(xmax-xmin)*rand(m,dimension)+xmin;%初始位置满足(-xmax,xmax)内均匀分布
        P=X;%P为每个粒子每代的最优位置
        last_f=f(X);
        [fmin,min_i]=min(last_f);%Pg为所有代中的最优位置 
        Pg=X(min_i,:);
        Savef(1)=fmin;
        N=0;
        for i=1:n
            v=w*v+c1*rand*(P-X)+c2*rand*(ones(m,1)*Pg-X);
            v=(v>vmax).*vmax+(v>=vmin & v<=vmax).*v+(v<vmin).*vmin;
            X=X+v*dt;
            X=(X>xmax).*xmax+(X>=xmin & X<=xmax).*X+(X<xmin).*xmin;
            new_f=f(X);%新的目标函数值
            update_j=find(new_f<last_f);
            P(update_j,:)=X(update_j,:);%修正每个粒子的历史最优值
            [new_fmin,min_i]=min(new_f);
            new_Pg=X(min_i,:);
            Pg=(new_fmin<fmin)*new_Pg+(new_fmin>=fmin)*Pg;
            last_f=new_f;%保存当前的函数值
            fmin=min(new_fmin,fmin);%更新函数最小值
             Savef(i)=fmin;
             if mod(i,floor(n/10))==0%10代记录一次种群参数
                 N=N+1;
                SaveData(:,:,N)=X;
             end
    %         w=w-i/n*0.7*w;
        end
        for j=1:10
            figure(j)
            plot(SaveData(:,1,j),SaveData(:,2,j),'o');
            xlim([-xmax,xmax])
            ylim([-xmax,xmax])
            axis tight
        end
        figure
        plot(Savef','b-')
        disp(Pg)
        disp(fmin)
    end
    

    这是一个粒子群算法的函数,可以求函数在一定范围内的极小值,下面我们来定义一个测试函数测试一下:
    f = ( 1 − x 1 ) 2 + ( x 1 2 − 2 x 2 ) 2 f=(1-x_1)^2+(x_1^2-2x_2)^2 f=(1x1)2+(x122x2)2
    我们很容易知道这个函数的极小值为 x ∗ = ( 1 , 0.5 ) T x^*=(1,0.5)^T x=(1,0.5)T
    我们将X和Y方向上的限制都设定为[-40,40],种群规模200,仿真代数100次

    f=@(x)(1-x(:,1)).^2+(x(:,1).^2-2*x(:,2)).^2;
    [Pg,fmin]=PSO(f,2,100,200,40,-40,40,-40)
    

    运行结果为:
    Alt

    序号粒子群位置随迭代次数的变化
    1Alt
    2在这里插入图片描述
    3在这里插入图片描述
    4在这里插入图片描述
    5在这里插入图片描述
    6在这里插入图片描述

    \qquad 读者们可以发现粒子群算法在迭代的过程中部分粒子甚至有反向运动的趋势,这是由于粒子的惯性速度过大导致的,在迭代点附近,惯性因子过大,会导致很长时间无法收敛到一个点,同时也很难找到相比上一次迭代更优的解。
    \qquad 我们只需要随着迭代次数增加减少我们的惯性因子 w w w即可,在这里作者给出的表达式为
    w k + 1 = w k − i / n ∗ 0.7 w k w_{k+1}=w_k-i/n*0.7w_k wk+1=wki/n0.7wk
    i i i为当前迭代代数, n n n为总代数。
    调用同样的测试函数,种群规模、迭代次数及其他参数不变,我们得到结果如下:
    在这里插入图片描述
    和上一次对比我们可以发现fmin小了很多,即在相同的计算规模下我们的最优解更加接近于极小值,可谓一个突破性的改变,我们再看看粒子随迭代次数的分布:

    序号粒子群位置随迭代次数的变化
    1在这里插入图片描述
    2在这里插入图片描述
    3在这里插入图片描述

    \qquad 可以发现粒子群的分布相比之前稳定得多,几乎全部收敛到了中心点,这便是减小惯性因子的作用。

    2.PID自整定

    2.1.基于M文件编写的PID参数自整定

    \qquad 利用M文件可以制作一个简单的PID引擎,用于测试我们的PSO算法是否有效,书写简单PID引擎的代码可以参考如下链接:
    【MATLAB】仿真PID控制
    \qquad 我们仅需要在此基础上做一些改变并做成并行算法即可,当然,我们还需要设置一个自适应度函数,在自适应度函数的作用下,PSO算法会选择到我们满意的参数。在调整PID参数时,我们通常会关注响应曲线的超调、上升时间、调节时间、峰值时间等等。利用这些参数组合为评价函数,让评价函数越小,仿真得到的参数越好,以此作为解的好坏的对比依据(相当于鸟离食物的距离)为了简便起见,我们只选取了调节时间 t s t_s ts和超调量 σ \sigma σ作为我们的比较依据,我们将评价函数设定为:
    f a c c e s s = l n ( t s 5 ⋅ 1 0 − 2 + 1 ) + l n ( σ 1 0 − 4 + 1 ) f_{access}=ln(\frac{t_s}{5\cdot 10^{-2}}+1)+ln(\frac{\sigma}{10^{-4}}+1) faccess=ln(5102ts+1)+ln(104σ+1)
    如果对该函数求导就可以发现:
    ∂ f ∂ t s ∣ t s = 5 ⋅ 1 0 − 2 = 1 2 ⋅ 5 ⋅ 1 0 − 2 \frac{\partial f}{\partial t_s}_{|t_s=5\cdot 10^{-2}}=\frac{1}{2}\cdot 5\cdot 10^{-2} tsfts=5102=215102
    ∂ f ∂ σ ∣ σ = 1 0 − 4 = 1 2 ⋅ 1 0 − 4 \frac{\partial f}{\partial \sigma}_{|\sigma=10^{-4}}=\frac{1}{2}\cdot 10^{-4} σfσ=104=21104
    \qquad 因此 f a c c e s s f_{access} faccess t s = 5 ⋅ 1 0 − 2 t_s=5\cdot 10^{-2} ts=5102的调整率为 2.5 × 1 0 − 2 2.5×10^{-2} 2.5×102,在 σ = 1 0 − 4 \sigma=10^{-4} σ=104时的调整率为 0.5 × 1 0 − 4 0.5×10^{-4} 0.5×104如果将期望参数定为 t s ∗ = 5 ⋅ 1 0 − 2 , σ ∗ = 1 0 − 4 t_s^*=5\cdot 10^{-2},\sigma ^*=10^{-4} ts=5102,σ=104则评价函数在小于任意一个参数的情况下就会下降,但是下降的速度会越来越大,而大于任意一个参数就会上升,但上升速率会越来越慢。由于两个参数均为正数,因此下降的速率不会大于 ∂ f ∂ t s ∣ t s = 0 = 1 \frac{\partial f}{\partial t_s}_{|t_s=0}=1 tsfts=0=1 ∂ f ∂ σ ∣ σ = 0 = 1 \frac{\partial f}{\partial \sigma}_{|\sigma=0}=1 σfσ=0=1,而上升的速率可以无限小。若将评价函数小视为自适应度高,则调整过后的系统参数总是在期望值附近调整,并尽量满足$所有参数都小于期望值的情况。因为在期望值附近上升速率是越来越慢而下降速率是越来越快的,而参数最小只能取到0。因此,有多个参数大于期望值而有少个参数小于期望值的和会比有少个参数大于期望值而多个参数小于期望值的和要大得多。
    \qquad 写好评价函数之后,我们可以设自适应度函数 F = − f a c c e s s F=-f_{access} F=faccess,因为我们的评价函数是越小越好的,而自适应度是越高越好。
    \qquad 根据1.2的算法步骤我们可以编写出粒子群算法,而PID并行仿真引擎的编写需要参考上文提供的连接。(simulink的PID仿真请参考下一小节2.2.)粒子群算法每迭代出一代的参数(PSO.m),就把这一代里面最好的参数 P G P_G PG返还给PID引擎进行仿真(PID_sim.m),再编写一个参数自识别的程序(parameter_cal.m),将需要识别的指标(调节时间、超调)计算出来,传递给 f a c c e s s f_{access} faccess进行计算(PID_access.m),将计算得到的参数

    【PID_sim.m】(PID并行仿真引擎)

    function [f_infty,tp,ts,sigma]=PID_sim(kp,ki,kd,debug)
        %kp,ki,kd为PID参数,T0为采样时间,total_t为仿真时间
        Tt=5e-4;%仿真采样周期
        T0=1e-2;%控制采样周期
        Tf=1e-3;%微分环节的滤波器系数
        alp=Tf/(Tt+Tf);
        total_t=1;
        N=floor(total_t/Tt);%样本总数
        M=numel(kp);%仿真序列数
        kp=reshape(kp,M,1);
        ki=reshape(ki,M,1);
        kd=reshape(kd,M,1);
        sys=tf(0.998,[0.021,1]);   %建立被控对象传递函数,即式4.1
        dsys=c2d(sys,Tt,'z');      %离散化
        [num,den]=tfdata(dsys,'v');   %
        e_1=zeros(M,1);      %前一时刻的偏差      
        Ee=zeros(M,1);       %累积偏差
        u_1=zeros(M,1);    %前一时刻的控制量
        y_1=zeros(M,1);       %前一时刻的输出
        ud_1=zeros(M,1);     %前一时刻的微分输出
        %预分配内存
        r=1500;
        y=zeros(M,N);
        u=zeros(M,N);
        for k=1:1:N
            y(:,k)=-1*den(2).*y_1+num(2)*u_1+num(1).*u(:,k);%系统响应输出序列
            e0=r-y(:,k);   %误差信号
            ud=alp.*ud_1+(1-alp)/T0.*kd.*(e0-e_1);
            u(:,k)=kp.*e0+T0*ki.*Ee+ud; %系统PID控制器输出序列
            Ee=Ee+e0;    %误差的累加和
            u_1=u(:,k);    	%前一个的控制器输出值
            y_1=y(:,k);    	%前一个的系统响应输出值
            e_1=e0;		%前一个误差信号的值
            ud_1=ud;
        end
        if debug %非调试模式下不显示也不打印图像
            plot(linspace(0,total_t,N),y)
        end
        [f_infty,tp,ts,sigma]=parameter_cal(y,Tt,2e-2,debug);%取得阶跃响应指标
    end
    

    【parameter_cal.m】(计算仿真曲线参数)

    function [f_infty,tp,ts,sigma]=parameter_cal(y,T0,delt_err,debug)
    %y是系统输出序列
    %T0是采样时间,可以结合时间序列点序号算出实际时间
    %delt_err是调节时间的误差精度
        M=size(y,1);N=size(y,2);%M为运算维度,N为时间序列长度
        f_infty=y(:,N);%稳态值序列
        err=y-f_infty*ones(1,N);%通过稳态值计算误差序列
        ferr=fliplr(abs(err));%倒序并取绝对值
        [~,ts_i]=max(ferr>delt_err*f_infty,[],2);
        ts_i=N*ones(M,1)-ts_i;
        ts=ts_i*T0;%调节时间
        [fp,tp]=max(y,[],2);%峰值和函数峰值
        tp=tp*T0;
        tp(abs(fp-f_infty)<1e-5)=NaN; %过阻尼无超调,没有峰值时间
        sigma=(fp-f_infty)./f_infty;
        if debug && M==1 %非调试模式下不显示
            disp(['系统稳态值为',num2str(f_infty)])
            disp(['系统超调量为',num2str(sigma*100),'%'])
            if isnan(tp)
                disp('系统不存在峰值时间')
            else
                disp(['系统峰值时间为',num2str(tp),'s'])
            end
            disp(['系统的调节时间为',num2str(ts),'s'])
        end
    end
    

    【PID_access.m】(根据参数得出评价度的函数,程序中为评价度越低越好)

    function y=PID_access(para)
    %PID调节性能的指标参数
    kp=para(:,1);ki=para(:,2);kd=para(:,3);
    [~,~,ts,sigma]=PID_sim(kp,ki,kd,false);
    y=log(ts/5e-2+1)+log(sigma/1e-2+1);
    end
    

    【PSO.m】(PSO算法优化函数)

    function [Pg,fmin]=PSO(dimension,n,m)
    %全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
        w=1;c1=2;c2=2;%速度惯性系数
        sigma_data=zeros(1,n);
        ts_data=zeros(1,n);
        dt=0.3;%位移仿真间隔
        vmax=1;%速度限幅
        xmax=[15,50,2];%位置限幅
        xmin=[0.2,0,0];
        v=zeros(m,dimension);%初始速度为0
        X=(xmax-xmin).*rand(m,dimension)+xmin;%初始位置满足(xmin,xmax)内均匀分布
        P=X;%P为每个粒子每代的最优位置
        last_f=PID_access(X);
        [fmin,min_i]=min(last_f);%Pg为所有代中的最优位置 
        Pg=X(min_i,:);
        N=0;
        for i=1:n
            v=w*v+c1*rand*(P-X)+c2*rand*(ones(m,1)*Pg-X);
            v=(v>vmax).*vmax+(v>=-vmax & v<=vmax).*v+(v<-vmax).*(-vmax);
            X=X+v*dt;
            X=(X>xmax).*xmax+(X>=xmin & X<=xmax).*X+(X<xmin).*(xmin);
            new_f=PID_access(X);%新的目标函数值
            update_j=find(new_f<last_f);
            P(update_j,:)=X(update_j,:);%修正每个粒子的历史最优值
            [new_fmin,min_i]=min(new_f);
            new_Pg=X(min_i,:);
            Pg=(new_fmin<fmin)*new_Pg+(new_fmin>=fmin)*Pg;
            last_f=new_f;%保存当前的函数值
            fmin=min(new_fmin,fmin);%更新函数最小值
            [~,~,ts,sigma]=PID_sim(Pg(1),Pg(2),Pg(3),true)
            sigma_data(1,i)=sigma;
            ts_data(1,i)=ts;
            hold on
        end
        legend('1','2','3','4','5','6','7','8','9','10','11','12','13','14','15')
        title('响应曲线随迭代次数的变化')
        figure
        plot(ts_data)
        title('调节时间随迭代次数的变化')
        figure
        plot(sigma_data)
        title('超调量随迭代次数的变化')
    end
    

    【main.m】(主函数,从这里开始)

    dimension=3;n=20;m=200;%PID为3位参数,n是迭代次数,m为种群规模
    [Pg,fmin]=PSO(dimension,n,m)
    

    仿真结果如下:
    在这里插入图片描述
    注:响应曲线可能要在左上角局部放大才能获得该图像。
    在这里插入图片描述
    在这里插入图片描述
    \qquad 这是一个简单的一阶惯性系统的PID调参仿真,用于模拟直流电机的转速模型,可以看出直流电机的转速响应是越来越快的(调节时间越来越小),超调量在第一次就达到了0,因此系统不存在超调。

    *2.2.复杂系统的PID自整定(基于simulink仿真)

    \qquad 我们以双闭环直流调速系统为例来说明PSO对PID自整定参数的作用。直流调速系统由外环(转速环)和内环(电流环)构成,转速环是一个PI调节器,调节转速为期望转速,用于抵抗负载扰动、电机参数变化带来的扰动;电流环也是一个PI调节器,与转速环不同的是,电流调节器的目的是在于让电流快速跟随给定值,以获得更好的动态特性。我们以电机的启动过程为例仿真这个系统,直流电机启动时,电流首先达到容许最大值,电机以最大加速度进行加速,到达期望转速时让电流迅速减小并达到负载转矩所需电流,电机维持期望转速不变。

    2.2.1.PSO优化PID的过程详解

    \qquad 下面我们设定我们需要达成的调节指标——由于启动时的电机电流会迅速上升到过载电流(允许的最大值),是一个恒定值,故会有一个电流平台,记上升到电流平台产生的超调为电流超调 σ i \sigma_i σi,从超调回到电流平台所需的时间为电流的调节时间 t s i t_{si} tsi,最终转速稳定前的超调记为转速超调 σ n \sigma_n σn,转速到达稳定值区域的时间为调节时间 t s n t_{sn} tsn
    【搭建simulink模型】
    双闭环调速系统的结构框图如下:
    在这里插入图片描述
    搭建simulink仿真模型如下:
    在这里插入图片描述
    注:上图可能不够清楚,根据以下参数搭建simulink模型,并且命名为double_circle.slx(和PSO的所有程序同名)
    电机的参数如下:

    • 额定电压/电流/功率:800V/2164A 1500KW
    • 额定转速:44r/min
    • 过载系数: λ i = 1.5 \lambda_i=1.5 λi=1.5
    • 励磁:110V/191A
    • 电枢电阻 R a R_a Ra:0.1 Ω \Omega Ω
    • 主回路电阻R:0.1014 Ω \Omega Ω
    • 电枢(主回路)电感L:3.25mH
    • 变流器放大倍数 K s K_s Ks=80
    • 变流器通态压降:3V
    • 系统飞轮惯量: 8.6 × 1 0 5 N ⋅ m 2 8.6×10^5N\cdot m^2 8.6×105Nm2

    \qquad 根据工程设计方法当然可以得到一套电流和转速的PI参数,计算的相关参数为: [ K n p , K n I , K i p , K i I ] E N = [ 7.3660 , 1.5658 , 0.8257 , 30.2773 ] [K_{np},K_{nI},K_{ip},K_{iI}]_{EN}=[7.3660 ,1.5658 ,0.8257, 30.2773] [Knp,KnI,Kip,KiI]EN=[7.3660,1.5658,0.8257,30.2773]
    通过仿真,我们得到的转速和电流曲线如下:

    额定负载启动空载启动
    在这里插入图片描述在这里插入图片描述

    调用我们的参数计算程序计算一下参数可以得到: σ n = 4.07 × 1 0 − 4 % , t s n = 0.7875 , σ i = 20.38 % , t s i = 0.0172 \sigma_n=4.07×10^{-4}\%,t_{sn}=0.7875,\sigma_i=20.38\%,t_{si}=0.0172 σn=4.07×104%,tsn=0.7875,σi=20.38%,tsi=0.0172
    \qquad 这个参数当然是比较理想的参数,如果我们不知道工程设计方法或者说工程设计方法不适用于我们的模型时,那么就可以使用PSO整定PID参数的方法。
    【修改simulink模型以适应PSO的M文件调用】
    \qquad 这里需要调整的参数为转速调节器ASR的 K n p , K n I K_{np},K_{nI} KnpKnI和电流调节器ACR的 K i p , K i I K_{ip},K_{iI} Kip,KiI,如果我们希望能实时调整simulink的参数进行动态仿真,我们需要直接将图中的几个比例环节依次设为符号( K n p , K n I , K i p , K i I K_{np},K_{nI},K_{ip},K_{iI} Knp,KnI,Kip,KiI),并在工作区实时修改这4个参数再进行仿真。如果我们要将实时调整PI参数的代码写在M文件中仿真,必须要将变量声明为全局(global),此时变量会在工作区里面,simulink就可以调用了。另存为PSO_double_circle.slx
    【编写辅助的M文件】
    cal_n.m:单组PID参数仿真结果的参数计算

    function [n_sigma,ts_n,I_sigma,ts_transient_i,n_err]=cal_n(parameters)
        global Kpn KIn Kpi KIi
        %将simulink所需参数导入工作区,以便仿真使用
        obj_slx='PSO_double_circle.slx';
        Kpn=parameters(1);
        KIn=parameters(2);
        Kpi=parameters(3);
        KIi=parameters(4);
        sim(obj_slx);
        I_data=I.signals.values;n_data=n.signals.values;time=I.time;
        n_infty=mean(n_data(numel(n_data)-10:numel(n_data)));
        n_err=abs(n_infty-10/0.227);%期望与输出的转速误差
        [i_max,tr_index]=max(I_data);
        diffI=[0;diff(I_data)];%电流微分
        pos=find(abs(diffI)<1);%找到电流波动较小的部分
        pos=pos(pos>tr_index);%排除启动时电流不变的情况
        min_pos=min(pos);
        I_transient=mean(I_data(min_pos+10:min_pos+20));%恒流升速阶段的暂稳态电流值
        I_sigma=(i_max-I_transient)/I_transient;
        ts_transient_i=time(min_pos-tr_index);%电流暂稳态的调节时间
        n_max=max(n_data);
        n_sigma=(n_max-n_infty)/n_infty;
        ts_n_index=numel(n_data)-find(abs(flip(n_data)-n_infty)/n_infty>2e-2,1);
        ts_n=time(ts_n_index);
    end
    

    PID_n_access.m:PID参数的仿真,以代为单位的

    function [y,ts_n,n_sigma,ts_transient_i,I_sigma]=PID_n_access(parameters_list)
        M=size(parameters_list,1);%计算需要仿真的组数(simulink无法并行仿真)
        ts_n0=0.8;n_sigma0=1e-2;I_sigma0=0.02;ts_i0=1e-4;%期望值
        %将simulink所需参数导入工作区,以便仿真使用
        y=zeros(M,1);
        n_sigma=zeros(M,1);ts_n=zeros(M,1);
        for sim_i=1:M
            parameters=parameters_list(sim_i,:);
            [n_sigma,ts_n,I_sigma,ts_transient_i,n_err]=cal_n(parameters);
            y(sim_i)=log(ts_n./ts_n0+1)+log(n_sigma./n_sigma0+1)+log(I_sigma/I_sigma0+1)+log(ts_transient_i/ts_i0+1)+exp(n_err*10);
        end
    end
    

    para_cal.m:参数计算,以代为单位的

    obj_slx='double_circle.slx';
    sim(obj_slx);
    I_data=I.signals.values;n_data=n.signals.values;time=I.time;
    I_infty=mean(I_data(numel(I_data)-10:numel(I_data)));
    n_infty=mean(n_data(numel(n_data)-10:numel(n_data)));
    [i_max,tr_index]=max(I_data);
    diffI=[0;diff(I_data)];%电流微分
    pos=find(abs(diffI)<1);%找到电流波动较小的部分
    pos=pos(pos>tr_index);%排除启动时电流不变的情况
    min_pos=min(pos);
    I_transient=mean(I_data(min_pos+10:min_pos+20));%恒流升速阶段的暂稳态电流值
    I_sigma=(i_max-I_transient)/I_transient;
    if I_sigma<1e-4
        tr_i=NaN;    
    else
        tr_i=time(tr_index);%电流的上升时间
    end
    ts_transient_i=time(min_pos-tr_index);%电流暂稳态的调节时间
    ts_i_index=numel(I_data)-find(abs(flip(I_data)-I_infty)/I_infty>2e-2,1);
    ts_i=time(ts_i_index);%转速稳态时电流的调节时间
    [n_max,tr_index]=max(n_data);
    n_sigma=(n_max-n_infty)/n_infty;
    %参数展示
    %转速参数展示
    if n_sigma<1e-4
        tr_n=NaN;
    else
        tr_n=time(tr_index);
    end
    ts_n_index=numel(n_data)-find(abs(flip(n_data)-n_infty)/n_infty>2e-2,1);
    ts_n=time(ts_n_index);
    if isnan(tr_n)
        disp('转速无超调')
    else
        disp(['转速上升时间为',num2str(tr_n),'s'])
        disp(['转速超调量为',num2str(100*n_sigma),'%'])
    end
    disp(['转速调节时间为',num2str(ts_n),'s'])
    %电流参数展示
    if isnan(tr_i)
        disp('电流无超调')
    else
        disp(['电流上升时间为',num2str(tr_i),'s'])
        disp(['电流超调量为',num2str(100*I_sigma),'%'])
        disp(['电流暂稳态调节时间为',num2str(ts_transient_i),'s'])
    end
    disp(['电流的最终调节时间为',num2str(ts_i),'s'])
    

    draw.m:仅根据一组PID参数画出转速波形图

    function draw(parameters)
        global Kpn Kin KpI KiI
        Kpn=parameters(1);
        Kin=parameters(2);
        KpI=parameters(3);
        KiI=parameters(4);
        sim('PSO_double_circle')
        n_data=n.signals.values;time=n.time;%获取转速序列及对应时刻
        plot(time,n_data)
    

    draw_n_I.m:根据一组PID参数绘制电流和转速曲线

    function draw_n_I(parameters)
        global Kpn KIn Kpi KIi
        Kpn=parameters(1);
        KIn=parameters(2);
        Kpi=parameters(3);
        KIi=parameters(4);
        sim('PSO_double_circle')
        n_data=n.signals.values;time=n.time;%获取转速序列及对应时刻
        I_data=I.signals.values;
        figure
        subplot(211)
        plot(time,I_data,'r-','LineWidth',1.5)
        title('电流波形')
        subplot(212)
        plot(time,n_data,'b-','LineWidth',1.5)
        title('转速波形')
    end
    

    PSO.m:PSO优化程序

    function [Pg,fmin]=PSO(n,m,xmax,xmin,cal_f,f)
    %全局粒子群算法,f为目标函数,dimension为维度,n为代数,m为种群规模
        dimension=numel(xmax);
        w=1;c1=2;c2=2;%速度惯性系数
        sigma_data=zeros(1,n);
        ts_data=zeros(1,n);
        dt=0.3;%位移仿真间隔
        vmax=0.1*(xmax-xmin);%速度限幅
        v=zeros(m,dimension);%初始速度为0
        X=(xmax-xmin).*rand(m,dimension)+xmin;%初始位置满足(xmin,xmax)内均匀分布
        P=X;%P为每个粒子每代的最优位置
        last_f=f(X);
        [fmin,min_i]=min(last_f);%Pg为所有代中的最优位置 
        Pg=X(min_i,:);
        last_Pg=Pg;
        legend_str=cell(0);
        legend_i=1;
        figure(1)
        legend_str{legend_i}=num2str(1);
        draw(Pg)
        for i=1:n
            v=w*v+c1*rand*(P-X)+c2*rand*(ones(m,1)*Pg-X);
            v=(v>vmax).*vmax+(v>=-vmax & v<=vmax).*v+(v<-vmax).*(-vmax);
            X=X+v*dt;
            X=(X>xmax).*xmax+(X>=xmin & X<=xmax).*X+(X<xmin).*(xmin);
            new_f=f(X);%新的目标函数值
            update_j=find(new_f<last_f);
            P(update_j,:)=X(update_j,:);%修正每个粒子的历史最优值
            [new_fmin,min_i]=min(new_f);
            new_Pg=X(min_i,:);
            Pg=(new_fmin<fmin)*new_Pg+(new_fmin>=fmin)*Pg;
            last_f=new_f;%保存当前的函数值
            fmin=min(new_fmin,fmin);%更新函数最小值
            [sigma,ts,~,~]=cal_f(Pg);
            sigma_data(1,i)=sigma;
            ts_data(1,i)=ts;
            if last_Pg~=Pg
                legend_i=legend_i+1;
                figure(1)
                legend_str{legend_i}=num2str(i);
                draw(Pg)
                hold on
            end
            last_Pg=Pg;
            disp(['迭代次数:',num2str(i)])
        end
        figure(1)
        legend(legend_str)
        title('响应曲线随迭代次数的变化')
        figure(2)
        plot(ts_data)
        title('调节时间随迭代次数的变化')
        figure(3)
        plot(sigma_data)
        title('超调量随迭代次数的变化')
    end
    

    set_para.m(设置参数程序,从这里开始)

    clear
    best_parameters=[8.4725,97,72,1.794,56.06];
    global Kpn KIn Kpi KIi
    Kpn=0;KIn=0;Kpi=0;KIi=0;
    n=10;m=100;xmax=[20,120,10,60];xmin=[0,0,0,0];cal_f=@cal_n;f=@PID_n_access;
    [Pg,fmin]=PSO(n,m,xmax,xmin,cal_f,f);
    
    

    我们将期望的参数设定为
    σ n ∗ = 0 , 1 % , t s n ∗ = 1 , σ i ∗ = 20 % , t s i ∗ = 0.005 \sigma_n^*=0,1\%,t_{sn}^*=1,\sigma_i^*=20\%,t_{si}^*=0.005 σn=0,1%,tsn=1,σi=20%,tsi=0.005利用PSO算法构造的评价函数为
    f a c c e s s = l n ( t s n 1 + 1 ) + l n ( σ n 0.001 + 1 ) + l n ( σ i 0.2 + 1 ) + l n ( t s i 0.005 + 1 ) ; f_{access}=ln(\frac{ts_n}{1}+1)+ln(\frac{\sigma_n}{0.001}+1)+ln(\frac{\sigma_i}{0.2}+1)+ln(\frac{ts_i}{0.005}+1); faccess=ln(1tsn+1)+ln(0.001σn+1)+ln(0.2σi+1)+ln(0.005tsi+1);
    得到的最佳参数为:
    [ K n p , K n I , K i p , K i I ] P S O = [ 18.2613 , 73.0965 , 1.6843 , 16.5703 ] P S O [K_{np},K_{nI},K_{ip},K_{iI}]_{PSO}=[18.2613,73.0965 ,1.6843,16.5703]_{PSO} [Knp,KnI,Kip,KiI]PSO=[18.2613,73.0965,1.6843,16.5703]PSO
    下的仿真曲线如下:
    超调量和调节时间随迭代次数变化
    在这里插入图片描述在这里插入图片描述
    PSO与工程设计方法电机启动曲线对比(额定负载启动)

    工程设计方法PSO方法
    在这里插入图片描述在这里插入图片描述

    通过计算程序得到的参数为: σ n = 6.7556 × 1 0 − 4 % , t s n = 0.7862 , σ i = 18.83 % , t s i = 0.01341 \sigma_n=6.7556×10^{-4}\%,t_{sn}=0.7862,\sigma_i=18.83\%,t_{si}=0.01341 σn=6.7556×104%,tsn=0.7862,σi=18.83%,tsi=0.01341
    工程设计方法的参数为:
    σ n = 4.07 × 1 0 − 4 % , t s n = 0.7875 , σ i = 20.38 % , t s i = 0.0172 \sigma_n=4.07×10^{-4}\%,t_{sn}=0.7875,\sigma_i=20.38\%,t_{si}=0.0172 σn=4.07×104%,tsn=0.7875,σi=20.38%,tsi=0.0172
    \qquad 可以发现除了最后一项指标之外,其余指标均达到了期望值要求,并且和工程设计方法得出的指标相比甚至还略优于它。

    2.2.2.在PSO优化过程中修改参数价值权重

    \qquad 转速的超调和调节时间,电流的超调和调节时间,我们最后评估的指标除了转速无静差的基本要求外,这4项指标就是我们的评价标准,如果我们更希望牺牲某一个或多个参数去优化另一个或另一些参数时,该怎么调节我们的算法呢?
    \qquad 答案也非常简单,调节的方法有2个,一个是期望值,一个是评价函数的各项系数。这里我们用调节期望值的方法做一个示范:
    \qquad 之前的期望值为:
    σ n ∗ = 0.1 % , t s n ∗ = 1 , σ i ∗ = 20 % , t s i ∗ = 0.005 \sigma_n^*=0.1\%,t_{sn}^*=1,\sigma_i^*=20\%,t_{si}^*=0.005 σn=0.1%,tsn=1,σi=20%,tsi=0.005
    修改期望值为:
    σ n ∗ = 10 % , t s n ∗ = 0.8 , σ i ∗ = 2 % , t s i ∗ = 0.0001 \sigma_n^*=10\%,t_{sn}^*=0.8,\sigma_i^*=2\%,t_{si}^*=0.0001 σn=10%,tsn=0.8,σi=2%,tsi=0.0001
    重新仿真,对比二者的结果如下:

    PSO1PSO2
    在这里插入图片描述在这里插入图片描述

    \qquad 可以看出转速超调量略微增大了一些,但是转速和电流的调节时间明显缩短,电流的超调量基本不变。这一点也表明多项优化指标是相互耦合的,可以牺牲一者满足另一者,只要参数在这方面是可调的。
    \qquad 需要注意的是,simulink中无法使用串行的算法,因此代码中使用时是串行的,可能仿真比较慢,请大家耐心等待(每迭代一次都会打印一次)
    \qquad 希望本文对您有帮助,谢谢阅读!

    展开全文
  • Simulink三相异步电机仿真

    万次阅读 多人点赞 2019-03-01 16:31:22
    前些时,利用CADe_SIMu V1.0仿真了一些常见常用的电机控制电路,但是CADe_SIMu V1.0软件只能仿真出电路的控制效果,而不能体现出元器件的相关参数、电路中电压、电流量的变化过程等。因此,本小学生利用这一周时间,...

    MATLAB/Simulink三相异步电机直接启动仿真
    小树不修不直溜,人不学习哏揪揪!隔了好长一段时间,本小学生终于又回归正常的学习与记录生活。前些时,利用CADe_SIMu V1.0仿真了一些常见常用的电机控制电路,但是CADe_SIMu V1.0软件只能仿真出电路的控制效果,而不能体现出元器件的相关参数、电路中电压、电流量的变化过程等。因此,本小学生利用这一周时间,并没有瞎扯蛋,而是学习了如何使用MATLAB中非常好用的电气仿真工具Simulink仿真最基本的电机启动环节,特总结此文,以留纪念,望大神指点!
    三相电机的直接启动
    三相电机的直接启动,也称作全压启动,是三相电机启动最简单的方法。启动时通过接触器将电机直接接入电网,设备简单,启动速度快且启动转矩大。直接起动适用于小容量电机带轻载的情况,对于大容量电机而言,这种起动方式具有十分显著的缺点,及起动电流较大,可达额定电流的4-7倍,易对电机和电网造成不利影响。因此,当额定功率P小于7.5kW时,可以采用直接起动。之前CADe_SIMu V1.0软件仿真的电机控制电路基本上都是选用的直接启动方式。
    Simulink三相异步电机直接启动仿真程序
    本小学生的电脑是windows7 32位系统,安装的MATLAB版本是2014a版本,这一版本的simulink与很多旧版本的simulink 大有不同,许多好用的仿真模块找不到了,而且有些模块的使用方法也有所不同,本校学生将以Simulink三相异步电机直接启动仿真程序为例,将创建仿真模型与数据参数修改的操作方法总结如下,方便大家参考。
    (1)创建Simulink仿真文件,点击“新建”下拉菜单中的“SIMULINK”再点击其下第1个选项“Simulink Model”,操作如图1所示,新建程序界面如图2所示。
    图1 Simulink程序创建操作
    图1 Simulink程序创建操作
    图2 Simulink程序创建界面
    图2 Simulink程序创建界面
    (2)元件的选择与放置:三相电源元件选择上,可供选择的模型有很多种,网上许多旧版本的simulink 仿真中直接利用3个AC Voltage Source模块构建三相电源,本小学生也对此种方案进行了尝试,不过AC Voltage Source模块搜索出的结果有两种,要选择如图3所示的右侧红色方框标注的器件,才能正常连接电机,本小学生也是在书写此文的过程中尝试实验成功的,下一篇文章中会体现出具体的仿真设计。
    图3 AC Voltage模块选择
    图3 AC Voltage模块选择
    本文中,小学生选用Three-Phase Programmable Voltage Source模块,作为仿真信号源,首先点击如图3所示的工具栏中图标,打开Simulink Library Browser,再在搜索栏中输入关键字Three-Phase,点击搜索栏右侧望远镜图标,搜索结果与Three-Phase Programmable Voltage Source模块位置如图4所示。
    图4 三相电源模块
    图4 三相电源模块
    三相电源模块的参数配置如图5所示,其中Amplitude(Vrms Ph-Ph)表示幅值,而Vrms表示有效值的含义,Phase表示相位角,Freq表示频率
    图5 三相电源参数配置
    图5 三相电源参数配置
    (3)三相电机选取:在搜索栏中输入关键字Machine,点击搜索栏右侧望远镜图标,搜索结果如图6所示,点击选择Asynchronous Machine如图中红色方框所示
    图6 三相异步电机模块
    图6 三相异步电机模块
    双击拖拽到面板上的三相电机模块,打开设置电机参数,设置Mechanical input为Torque TM(转矩),Rotor Type为Squirrel-cage(鼠笼式)如图7所示。电机的具体电气参数设置如图8所示,第一行中Nominal power表示电机的额定功率,其他两个分别是供电幅值和频率,第二行Stator是电机定子,而resistance和inductance相信许多学过电子类知识的小伙伴,都能通过参数的单位猜出,是电阻与电感参数设置。第三行Rotor是电机转子的参数设置,我设置的电机参数如图8所示。
    图7 电机参数设置1
    图7 电机参数设置(1)
    图8 电机参数设置2
    图8 电机参数设置(2)
    由于我在设置电机机械输入量时,设置了转矩量,因此需要如图9所示的Constant(连续量)模块,其参数设置如图10所示,注意图10中红色方框标注的选项对勾一定要去掉,不然输出的会是向量,而不是常量。
    图9 连续量模块
    图9 连续量模块
    图10 连续量参数配置
    图10 连续量参数配置
    (4)连线操作与CADeV1.0等其他仿真软件基本相似,都是选好连线的起始位置后,左键按住拖拽到连线终端,松开按键即可。所需注意的是,当连线显示是红色时表示连线失败,而且必须从起点开始至终点结束,不能逐步连接如图11所示即为失败情况,图12 所示是连线成功的效果。
    图11 接线失败情况
    图11 接线失败情况
    图12  接线成功情况
    图12 接线成功情况
    (5)检测模块:MATLAB2014软件版本的simulink 中没有了Machines Measurement Demux模块。所以为了检测电机的一些运行参数,本小学生选用总线Bus Selector模块,将电机的参数分离出来,如图13所示选取Bus Selector模块。双击模块,选择测量的电机参数如图14,我这里检测转子、定子电流、电机转速、转矩四组检测量如图15 所示。
    图13 总线模块选择
    图13 总线模块选择
    图14 总线参数选择
    图14 总线参数选择
    图15 检测变量选择
    图15 检测变量选择
    配置好参数后的Bus Selector会如图16所示多出2个支路出来
    图16 配置好的总线模块
    图16 配置好的总线模块
    (6)最后配置示波器模块,搜索关键字Scope,如图17所示为选择操作过程。双击Scope模块,调整示波器的一些参数。如图18所示,点击齿轮图标,进入参数配置界面General下,配置Number of axes(坐标数目)与被测参量数目相同,设为4。TIME range 为横坐标最大值,我这里设置成1,Tick labels表示坐标标签,实际上是选择各个波形图是否显示横坐标值,我这里选择all,如图19所示。第2个界面History是可以选择是否存储测量数据,第三个界面Style可以设置波形图颜色,这里将背景设置成白色线设置成黑色,如图20所示。
    图17 示波器模块
    图17 示波器模块
    图18 示波器参数设置1
    图18 示波器参数设置(1)
    图19 示波器参数设置2
    图19 示波器参数设置(2)
    调试好参数的示波器效果如图20所示
    图20 示波器调试效果
    图20 示波器调试效果
    (7)配置POWERGUI模块,POWERGUI是电力系统仿真模块,如果你是学习电气工程等电类专业则可以用到它,可进行电网稳定性分析、傅里叶分解、潮流计算、阻抗频率响应等计算。选取操作如图21所示。虽然并没有用上,但是缺少此模块仿真会失败,不明原因,有了解原因的小伙伴或大神可评论中赐教,谢谢!
    图21 POWERGUI模块
    图21 POWERGUI模块
    (8)最后连接好各个模块如图22所示
    图22 完整的仿真程序
    图22 完整的仿真程序
    仿真文件保存与仿真结果展示
    (1)保存操作:一定要保存,不然就全白做了。首先点击File下拉菜单选择Save,或CTRL+S如图23所示,弹出save as界面如图24所示选择保存路径,保存文件名一定是英文名,不然保存不成功提示如图25所示。
    图27 保存操作
    图27 保存操作
    图28 保存正确设置
    图28 保存正确设置
    图28 保存错误提醒
    图28 保存错误提醒
    (2)仿真操作如图29所示点击图中红色方框的开始按钮,即可开始仿真,本程序仿真速度特别快。仿真完成后,双击Scope模块,展示出我们期望的参量波形,起始时电流会有短暂的急剧变化,平稳后几乎不变,平稳后的稳定转速为1500r/min,如图30所示,即为仿真测量的四组数据对应波形。
    图29 仿真操作
    图29 仿真操作
    图30 仿真数据波形
    图30 仿真数据波形
    最后分享此工程的仿真文件和图片给各位老铁:
    百度网盘永久有效链接:
    链接:https://pan.baidu.com/s/1EwaIiMlM3j1H-clK33Y-WQ
    提取码:9rgc
    良心博主,且看且珍惜,如需分享,表明转载,小弟不胜感激!

    展开全文
  • 四足机器人运动控制第一章 序第二章 运动状态姿态控制运动控制第三章 步态第四章 CPG震荡单元模型介绍Cpg模型分类基于HOPF振荡器的Cpg单元模型参考文献 第一章 序 足式机器人较传统的四轮式和履带式有着无与伦比的...
  • 第五章 卷积神经网络

    千次阅读 多人点赞 2020-01-13 07:49:37
    第五章 卷积神经网络第五章 卷积神经网络卷积一维卷积二维卷积互相关卷积的变种卷积的数学性质交换性导数卷积神经网络用卷积来代替全连接卷积层汇聚层(池化层)典型的卷积网络结构参数学习误差项的计算几种典型的卷积...
  • 拓端tecdat|R语言时变参数VAR随机模型

    千次阅读 2019-06-13 22:06:27
    时变参数VAR随机模型是一种新的计量经济学方法,用于在具有随机波动率和相关状态转移的时变参数向量自回归(VAR)的大模型空间中执行随机模型规范搜索(SMSS)。这是由于过度拟合的关注以及这些高度参数模型中通常...
  • 不知道小伙伴使用pytorch时是否遇到了下面的情况:我们加载训练好的参数,然后输入经过网络,但是不经过梯度下降训练,也就是说模型参数没有更改,这是我们仍然会发现我们输出的指标(例如准确率)下降。...
  • SPSS(五)SPSS之相关分析与线性回归模型(图文+数据集) 在讲解线性回归模型之前,先来学习相关分析的知识点,因为相关分析与回归有着密切的联系 相关分析 任意多个变量都可以考虑相关问题,不单单局限于两个...
  • 由于尺度不变性,该模型仅具有两个独立的参数,并且对于受DM残留物密度限制的参数空间,可以发生强烈的一阶电弱相变。 在此模型中,对于参数空间的一小部分,DM核子横截面低于中微子底限,因此,将来的直接检测实验...
  • Stata: 空间面板数据模型及Stata实现

    万次阅读 多人点赞 2019-05-10 10:37:56
    由于面板数据模型所具有的众多优点 (刻画个体异质性,减弱模型共线性和增加自由度等),其被广泛应用于实证计量中。在 「 Stata: 面板数据模型-一文读懂」 文中,我们已对面板数据模型进行了介绍。 然而,当研究样本...
  • 我们表明,用于描述暗物质与标准模型粒子之间相互作用的简化模型在一般情况下并不符合规范不变性,并且在参数空间的较大区域中可能会破坏微扰统一性。 解决这些矛盾所必需的修改可能意味着现象学更加丰富,并导致对...
  • 相机模型(内参数,外参数

    千次阅读 2018-04-26 21:05:54
    相机成像的过程实际是将真实的三维空间中的三维点映射到成像平面(二维空间)过程,可以简单的使用小孔成像模型来描述该过程,以了解成像过程中三维空间到二位图像空间的变换过程。 本文包含两部分内容,首先介绍...
  • 在含水量保持不变的情况下,邓肯一张E~μ模型的8个参数中,参数K、n、Rr、C、φ、G均随压实度的提高而增大;F、D随K的变化规律不明显。在含水量和压实度保持不变的情况下进行三轴试验,土样的(σ1-σ3)f、(σ1-...
  • 也可以划分为三个互斥的集合,此时就增加一个验证集,用于调试参数和选择模型 。 直接采用 sklearn 库的 train_test_split 函数即可实现,一个简单的示例代码如下,这里简单调用 knn 算法,采用 Iris 数据...
  • 去年TensorFlow官方推出了模型优化工具,最多能将模型尺寸减小4倍,运行速度提高3倍。 最近现又有一款新工具加入模型优化“豪华套餐”,这就是基于Keras的剪枝优化工具。 训练AI模型有时需要大量硬件资源,但不是...
  • 1 弱磁扩速理论 PMSM弱磁的思想来源于他励直流电动机的调磁... 永磁同步电机的励磁磁通是由永磁体提供的,这个磁通是恒定不变的。这个时候如果我们想降低磁通强度,就只能通过增大定子电流的去磁分量来削弱气隙...
  • 在以Hooft Rξ规量化的U(1)希格斯模型中研究了一组局部规范不变复合算符的光谱特性。 这些算子使我们能够对理论的频谱进行规范不变的描述,从而在使用标准基本场时超过了某些商品。 一环评估相应的两点相关函数,...
  • 机器学习中的模型参数估计方法:极大似然估计、贝叶斯估计、最大后验估计。
  • 深度学习 模型训练超参数调整总结

    千次阅读 2020-04-01 17:23:02
    深度学习 模型训练超参数调整总结 在深度神经网络中,超参数的调整是一项必备技能,通过观察在训练过程中的监测指标如损失loss和准确率来判断当前模型处于什么样的训练状态,及时调整超参数以更科学地训练模型能够...
  • 面板数据分析步骤及流程-R语言

    万次阅读 多人点赞 2016-08-16 16:49:55
    面板数据模型选择及分析步骤;附R语言代码
  • 作者用 YOLOv3 做人手检测(在 oxford hand 数据集上训练的),并进行了模型剪枝,剪枝后YOLOv3 模型参数量减少 80% ,FLOPs 降低 70%,推断的速度提高了100%,而 mAP 基本保持不变! 可谓是相当成功了! 剪枝前后...
  • 鲁棒性

    千次阅读 2018-05-23 00:34:31
    鲁棒性(robustness:健壮性|抗变换性):统计学中的专业术语,自20实际70年代初开始在控制理论的研究中流行起来,用来体现控制系统在一定(结构,大小)特性或参数扰动|摄动的情况下维持系统的某些满意性能的不敏感性|...
  • 最近在做一个用简单的lstm预测时序数据分类的模型,在训练过程中,出现了loss在持续减小,而acc和val_acc始终不变的情况,经过不断的调参,找到了其中的问题,在此记录下: 问题描述: 具体情况如下: Epoch ...
  • 模型采用简单的四层linear以及Relu、Sigmoid,实现二分类问题 loss采用的是交叉熵和Focal loss(测试Focal loss性能) 优化方式采用的是Adam+StepLR 二、LOSS不变的原因 1. 背景:训练集和测试集的loss都不变了、训练...
  • 以该输出不变集作为多参数规划问题中状态区域约束限制的初始条件,通过反复求解多参数规划问题和不断改变状态区域约束限制,能够逐步扩大显式模型预测控制系统的无限时间可行区域,直到可行域不再继续扩大.算法收敛时...
  • 自适应控制基本思想

    万次阅读 多人点赞 2017-09-22 10:18:20
    但实践中我们还会遇到结构和参数都未知的对象,比如一些运行机理特别复杂,目前尚未被人们充分理解的对象,不可能建立有效的数学模型,因而无法沿用基于数学模型的方法解决其控制问题,这时需要借助人工智能学科,也...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 132,987
精华内容 53,194
关键字:

不变参数模型