位置式\增量式PID、模糊PID、BRF-PID的Matlab实现以及封装

位置式\增量式PID、模糊PID、BRF-PID的Matlab实现以及封装

简要

PID控制由于结构简单，控制效果好，在工业上有广泛的应用，从出现到至今已经有一百多年的历史；PID控制本质是一种线性控制器，它是比例（P）、积分（I）、微分（D）的线性组合，所以一旦参数确定，控制效果就已经确定，在面对非线性系统时候，控制效果不是太理想。所谓的智能PID就是通过把智能算法与PID相结合，各取所长得到智能PID算法，其算法本质就是通过规则或者学习算法来在线的改变PID的三个参数，从而面对时变、非线性系统时，能有较强的鲁棒性。

本文主要包括位置式\增量式PID、模糊PID、BRF-PID的Matlab实现以及封装，其中包括了完整的代码。Matlab也是一门面向对象语言，利用面向对象的强大功能，本文通过策略模式把算法经行了封装，方便以后使用。以下是完整的代码实现：

Model类

此类负责储存数据

classdef Model < handle
% 策略模式：定义一系列的算法，把它们一个个封装起来，并且使它们可以互相替换
% 此类对象中存放数据，而在Strategy中存放算法，此类可以选择不同的算法
    properties
         sampleTime      % 采样时间
         simTime         % 仿真时间
         rin             % 输入
         sys_Plant       % 被控对象传递函数
         dsys            % 离散模型
         den             % 离散模型分母
         num             % 离散信号分子
         ui
         yi
         i
         strategyObj     % 策略对象
    end

    methods 
        function obj = Model(sampleTime,simTime,rin,plant)
            % 初始化数据。
            obj.initialize(sampleTime,simTime,rin)
            % 建立模型及初始化其他参数
            obj.buildPlant(plant);
        end
 
        function initialize(obj,sampleTime,simTime,rin)
            obj.sampleTime = sampleTime;
            obj.simTime = simTime;
            obj.i =  obj.simTime/obj.sampleTime;
            if strcmp(rin,'step')
                obj.rin = ones(obj.i,1);
            elseif strcmp(rin,'impulse')
                obj.rin = zeros(obj.i,1);
                obj.rin(1) = 1/obj.sampleTime;
            elseif strcmp(rin,'square wave')
                x = 1:obj.i;
                obj.rin = sign(sin(2*pi*obj.sampleTime.*x));
            end
        end
        
        % 被控对象传递函数建立及其他参数初始化。
        function buildPlant(obj,plant)
            s = tf('s');
            obj.sys_Plant = eval(plant);
            % 零阶保持器离散化
            obj.dsys = c2d(obj.sys_Plant,obj.sampleTime,'z');
            [obj.num,obj.den] = tfdata(obj.dsys,'v');
            % 其他初始化
            obj.initialize_other();
        end

        function initialize_other(obj)
            n = length(obj.den);
            obj.ui = zeros(n-1,1);
            obj.yi = zeros(n-1,1);
        end

        function progress(obj)
            obj.initialize_other()
            obj.strategyObj.progress(obj);
        end
    end
end

Strategy 类

此类为抽象类，定义算法接口

classdef Strategy < handle
% 策略模式
% 策略抽象基类，定义算法接口，有PID_Strategy,FuzzyPID_Strategy,BRF_PID_Strategy,SingleNeuralPID_Strategy等子类。
    properties
         PID_param = [0 0 0];   % 向量依次代表pid三个参数
         yout                   % 输出
         stepTime               % 离散时间
         stepRin                % 离散输入
         error                  % 误差变化
         PID                    % PID三个参数变化过程
    end
    methods
        function obj = Strategy(P,I,D)
            obj.PID_param = [P I D];
        end
    end
    methods(Abstract)
        progress(obj) 
    end
end

PID_Strategy类

此类为位置式、增量式PID算法的具体实现。

classdef PID_Strategy < Strategy
    properties
        displacementPID             % 位置式PID和增量式PID切换，0代表增量式PID，1代表位置式PID
    end
    methods
        function obj = PID_Strategy(P,I,D,displacementPID)
            obj = obj@Strategy(P,I,D);
            obj.displacementPID = displacementPID;
        end

        function progress(obj,modelObj)
            x = [0 0 0]';
            i = modelObj.i;
            obj.stepRin = modelObj.rin;
            obj.yout = zeros(i,1);
            time = zeros(i,1);
            % 位置式PID
            u = zeros(i,1);                 % 控制器输出
            e = zeros(i,1);                 % 误差
            e_1 = 0;                        % 上一时刻误差
            % 增量式PID
            du = u;
            u_1 = 0;
            e_2 = 0;
            for k = 1:i
                time(k) = k*modelObj.sampleTime;
                % 输出
                obj.yout(k) = modelObj.num(2:end)*modelObj.ui - modelObj.den(2:end)*modelObj.yi;
                % 误差
                e(k) = obj.stepRin(k) - obj.yout(k);
                % 位置式PID算法
                if obj.displacementPID 
                    u(k) = obj.PID_param*x;
                    if u(k) >= 10
                        u(k) = 10;
                    end
                    if u(k) <= -10
                        u(k) = -10;
                    end
                    x(1) = e(k);
                    x(2) = x(2) + e(k)*modelObj.sampleTime; 
                    x(3) = e(k) - e_1;
                % 增量式PID算法
                else
                    du(k) = obj.PID_param*x;
                    u(k) = u_1 + du(k);
                    if u(k) >= 10
                        u(k) = 10;
                    end
                    if u(k) <= -10
                        u(k) = -10;
                    end
                    x(1) = e(k) - e_1;
                    x(2) = e(k);
                    x(3) = e(k) - 2*e_1 + e_2; 
                    u_1 = u(k);
                    e_2 = e_1;
                end              
                % 跟新数据
                modelObj.ui(2:end) = modelObj.ui(1:end-1);
                modelObj.ui(1) = u(k);
                modelObj.yi(2:end) = modelObj.yi(1:end-1);
                modelObj.yi(1) = obj.yout(k); 
                e_1 = e(k);
            end
            obj.stepTime = time;
            figure
            plot(time,obj.yout)
            grid on
            obj.error = e;
        end
    end
end

Fuzzy_PID_Strategy类

此类为模糊PID的具体实现

classdef Fuzzy_PID_Strategy < Strategy
% 模糊PID控制算法。二输入，三输出
    properties
        fid 				            % 模糊逻辑控制文件
    end
    
    methods
        function obj = Fuzzy_PID_Strategy(P,I,D,file)
            obj = obj@Strategy(P,I,D);
            obj.fid = readfis(file);
        end

        function progress(obj,modelObj)
            x = [0 0 0]';
            i = modelObj.i;
            obj.PID = zeros(i,3);
            obj.stepRin = modelObj.rin;
            obj.yout = zeros(i,1);
            time = zeros(i,1);
            u = zeros(i,1);                 % 控制器输出
            e = zeros(i,1);                 % 误差
            e_1 = 0;                        % 上一时刻误差
            ec_1 = 0;                       % 上一时刻误差变化           
            for k = 1:i
                time(k) = k*modelObj.sampleTime;
                % 根据模糊规则输出PID变化
                deltaPID = evalfis(obj.fid,[e_1,ec_1]);
                obj.PID(k,:) = obj.PID_param + deltaPID;
                u(k) = obj.PID(k,:)*x;
                if u(k) >= 10
                    u(k) = 10;
                end
                if u(k) <= -10
                    u(k) = -10;
                end
                % 输出
                obj.yout(k) = modelObj.num(2:end)*modelObj.ui - modelObj.den(2:end)*modelObj.yi;
                % 误差
                e(k) = obj.stepRin(k) - obj.yout(k);
                x(1) = e(k);
                x(2) = x(2) + e(k)*modelObj.sampleTime; 
                x(3) = e(k) - e_1; 
                % 跟新数据
                modelObj.ui(2:end) = modelObj.ui(1:end-1);
                modelObj.ui(1) = u(k);
                modelObj.yi(2:end) = modelObj.yi(1:end-1);
                modelObj.yi(1) = obj.yout(k); 
                e_1 = e(k);
                ec_1 = x(3);
            end
            obj.stepTime = time;
            figure
            plot(time,obj.yout)
            grid on
            obj.error = e;
        end
    end
end

BRF_PID_Strategy类

此类为径向基函数PID的具体实现

classdef BRF_PID_Strategy < Strategy
% 径向基函数网络PID算法,核函数为高斯基函数，网络为3-6-1形式
    properties
        alpha                            % RBF网络学习效率
        eta                              % RBF网络动量因子
        b = 4                            % 高斯基函数宽度
        c                                % 高斯基函数中心矩阵
        hide_num = 6                     % 隐含层个数
        out_range                        % 模型输出范围
        etaPID = [0.2 0.2 0.2]           % PID三个参数学习效率
        ymout                            % BRF网络逼近输出
        dyout                            % BRF网络雅可比阵
        displacementPID                  % 位置式PID和增量式PID切换，0代表增量式PID，1代表位置式PID
    end
    
    methods
        function obj = BRF_PID_Strategy(P,I,D,etaPID,alpha,eta,b,out_range,hide_num,displacementPID)
            obj = obj@Strategy(P,I,D);  
            obj.etaPID = etaPID;
            obj.alpha = alpha;
            obj.eta = eta;
            obj.b = b;
            obj.out_range = out_range;
            obj.hide_num = hide_num;
            obj.displacementPID = displacementPID;
            obj.c = [linspace(out_range(1),out_range(2),obj.hide_num);
                    linspace(out_range(1),out_range(2),obj.hide_num);
                    linspace(out_range(1),out_range(2),obj.hide_num)];
        end

        function progress(obj,modelObj)
            x = [0 0 0]';  
            y_1 = 0;                        % 上一时刻输出
            % 初始化神经网络各参数
            w = zeros(obj.hide_num,1);      % 初始权矩阵
            w_1 = w;
            w_2 = w_1;
            h = zeros(obj.hide_num,1);      % 初始化高斯基函数输出
            RBF_input = [0 0 0]';           % 初始化神经网络输入
            % 初始化模型输入输出参数
            i = modelObj.i;
            obj.PID = zeros(i,3);
            obj.stepRin = modelObj.rin;
            obj.yout = zeros(i,1);
            obj.ymout = zeros(i,1);
            obj.dyout = zeros(i,1);
            time = zeros(i,1);
            % 位置式PID
            u = zeros(i,1);                 % 控制器输出
            e = zeros(i,1);                 % 误差
            e_1 = 0;                        % 上一时刻误差
            % 增量式PID
            du = u;
            u_1 = 0;
            e_2 = 0;
            
            for k = 1:i
                time(k) = k*modelObj.sampleTime;
                % 输出
                obj.yout(k) = modelObj.num(2:end)*modelObj.ui - modelObj.den(2:end)*modelObj.yi;
                for j = 1:obj.hide_num
                    h(j) = exp(-norm(RBF_input - obj.c(:,j))^2/(2*obj.b^2));
                end
                obj.ymout(k) =  w'*h;
                d_w = obj.alpha*(obj.yout(k) - obj.ymout(k))*h;
                w = w_1 + d_w + obj.eta*(w_1 - w_2);
                yu = w.*h.*(-x(1) + obj.c(1,:))'/obj.b^2;
                obj.dyout(k) = sum(yu);
                % 误差
                e(k) = obj.stepRin(k) - obj.yout(k);
                obj.PID_param = obj.PID_param + e(k)*obj.dyout(k)*x'.*obj.etaPID;
                obj.PID(k,:) = obj.PID_param;
                % 位置式PID算法
                if obj.displacementPID
                    u(k) = obj.PID_param*x;
                    if u(k) >= 10
                        u(k) = 10;
                    end
                    if u(k) <= -10
                        u(k) = -10;
                    end
                    x(1) = e(k);
                    x(2) = x(2) + e(k)*modelObj.sampleTime; 
                    x(3) = e(k) - e_1;
                    RBF_input = [u(k) obj.yout(k) y_1]';
                % 增量式PID算法
                else
                    du(k) = obj.PID_param*x;
                    u(k) = u_1 + du(k);
                    if u(k) >= 10
                        u(k) = 10;
                    end
                    if u(k) <= -10
                        u(k) = -10;
                    end
                    x(1) = e(k) - e_1;
                    x(2) = e(k);
                    x(3) = e(k) - 2*e_1 + e_2; 
                    RBF_input = [du(k) obj.yout(k) y_1]';
                    u_1 = u(k);
                    e_2 = e_1;
                end
                % 跟新数据
                modelObj.ui(2:end) = modelObj.ui(1:end-1);
                modelObj.ui(1) = u(k);
                modelObj.yi(2:end) = modelObj.yi(1:end-1);
                modelObj.yi(1) = obj.yout(k); 
                e_1 = e(k);
                w_2 = w_1;
                w_1 = w;
                y_1 = obj.yout(k);
            end
            obj.stepTime = time;
            obj.error = e;
            figure
            plot(time,obj.yout)
            grid on
            figure
            plot(time,obj.yout,'r',time,obj.ymout,'b');
        end
    end
end

测试

测试对象为三阶的一型系统，采样时间$T_S=0.001s$，仿真时长$T=0.5s$,输入为单位越阶。

%%  PID算法/增量式PID算法调试
a = Model(0.001,0.5,'step','523500/(s^3+87.35*s^2+10470*s)');
a.strategyObj = PID_Strategy(1,0,0,1);
% a.strategyObj = PID_Strategy(1,0,0,0);
a.progress()
%%  神经网络PID/神经网络增量式PID算法调试
% a = Model(0.001,0.5,'step','523500/(s^3+87.35*s^2+10470*s)');
% isa(a.strategyObj,'Strategy')
% a.strategyObj = BRF_PID_Strategy(0.8,0,0,[0.9,0.9,0.2],0.4,0.01,2,[-1,1],6,1);
% a.progress()