自适应模糊PID控制，理论+仿真

海雅信子

已于 2024-05-31 18:15:20 修改

阅读量1.1k

点赞数 9

分类专栏：控制文章标签：算法 matlab

于 2024-05-31 14:56:06 首次发布

本文链接：https://blog.csdn.net/adxmuying/article/details/139351987

版权

控制专栏收录该内容

4 篇文章 0 订阅

订阅专栏

不太喜欢写理论，模糊PID的相关内容可以参考书籍《模糊控制系统及应用》，《先进PID控制MATLAB仿真》里也简单描述了一下基本理论。
简单来说，自适应模糊PID控制器以误差 $e$ 和误差变化 $e c$ 作为模糊控制器的输入，利用模糊控制规则在线对PID参数进行调整,以满足不同时刻的 $e$ 和 $e c$ 对PID参数自整定的要求。自适应模糊PID控制器系统结构如下图所示：
在这里插入图片描述
控制率和PID一样：
$k_p e(k)+k_i T \sum^k e(j)+k_d \frac{1}{T}(e(k)-e(k-1))$
$\delta kp, \delta ki, \delta kd$ 模糊规则如下：

输入输出的隶属函数如下：

大致流程如下：

将 $e$ 和 $e c$ 值量化为论域值。
利用该论域值计算PID参数增量 $\delta kp, \delta ki, \delta kd$ 的论域值。
将PID参数增量论域值反量化为实际的增量值。
自调整PID参数。

假设控制对象为：
$\frac{133}{s^2+25\cdot s}$
代码如下（看懂就都会了）：
FuzzyPID.m

% created by hyacinth on 2024/5/31
classdef FuzzyPID < handle
    properties
        e_max       % 模糊集论域
        e_min
        ec_max
        ec_min
        delta_kpmax % pid参数增量最大最小值
        delta_kpmin
        delta_kimax
        delta_kimin
        delta_kdmax
        delta_kdmin
        dt

        kp          % pid参数
        ki
        kd

        qdelta_kp   % pid参数增量论域值
        qdelta_ki
        qdelta_kd

        delta_kp    % pid参数增量值
        delta_ki
        delta_kd

        qerro       % 误差相关论域值
        qerro_c

        pre_err     % 前一时刻误差
        erri        % 积分项

        threshold

        e_gradmembership        % 输入的隶属度,形状(1,2) or (1,1)
        ec_gradmembership
        e_grad_idx              % 隶属度在模糊子集中的位置,形状(1,2) or (1,1)
        ec_grad_idx
        delta_kp_grad_sums      % pid参数增量总隶属度,形状(1,7)
        delta_ki_grad_sums
        delta_kd_grad_sums

        kp_rule_list            % 模糊规则,形状(7,7)
        ki_rule_list
        kd_rule_list

        e_membership_values     % 输入输出的模糊子集,形状(1,7)
        ec_membership_values
        kp_menbership_values
        ki_menbership_values
        kd_menbership_values
    end
    methods
        % 构造函数，初始化模糊系统
        function obj = FuzzyPID(e_max,e_min,ec_max,ec_min,kp0,ki0,kd0, ...
                delta_kpmax,delta_kpmin,delta_kimax,delta_kimin,delta_kdmax,delta_kdmin,threshold,dt)
            obj.e_max = e_max;
            obj.e_min = e_min;
            obj.ec_max = ec_max;
            obj.ec_min = ec_min;
            obj.delta_kpmax = delta_kpmax;
            obj.delta_kpmin = delta_kpmin;
            obj.delta_kimax = delta_kimax;
            obj.delta_kimin = delta_kimin;
            obj.delta_kdmax = delta_kdmax;
            obj.delta_kdmin = delta_kdmin;
            
            obj.kp = kp0;
            obj.ki = ki0;
            obj.kd = kd0;
            obj.qdelta_kp = 0;
            obj.qdelta_ki = 0;
            obj.qdelta_kd = 0;
            obj.pre_err = 0;
            obj.erri = 0;
            obj.threshold = threshold;
            obj.dt = dt;

            NB = -3;
            NM = -2;
            NS = -1;
            ZO = 0;
            PS = 1;
            PM = 2;
            PB = 3;

            obj.e_membership_values = [NB,NM,NS,ZO,PS,PM,PB];
            obj.ec_membership_values = [NB,NM,NS,ZO,PS,PM,PB];
            obj.kp_menbership_values = [NB,NM,NS,ZO,PS,PM,PB];
            obj.ki_menbership_values = [NB,NM,NS,ZO,PS,PM,PB];
            obj.kd_menbership_values = [NB,NM,NS,ZO,PS,PM,PB];

            obj.kp_rule_list = [PB,PB,PM,PM,PS,ZO,ZO;
                                PB,PB,PM,PS,PS,ZO,NS;
                                PM,PM,PM,PS,ZO,NS,NS;
                                PM,PM,PS,ZO,NS,NM,NM;
                                PS,PS,ZO,NS,NS,NM,NM;
                                PS,ZO,NS,NM,NM,NM,NB;
                                ZO,ZO,NM,NM,NM,NB,NB];
            obj.ki_rule_list = [NB,NB,NM,NM,NS,ZO,ZO;
                                NB,NB,NM,NS,NS,ZO,ZO;
                                NB,NM,NS,NS,ZO,PS,PS;
                                NM,NM,NS,ZO,PS,PM,PM;
                                NM,NS,ZO,PS,PS,PM,PB;
                                ZO,ZO,PS,PS,PM,PB,PB;
                                ZO,ZO,PS,PM,PM,PB,PB];
            obj.kd_rule_list = [PS,NS,NB,NB,NB,NM,PS;
                                PS,NS,NB,NM,NM,NS,ZO;
                                ZO,NS,NM,NM,NS,NS,ZO;
                                ZO,NS,NS,NS,NS,NS,ZO;
                                ZO,ZO,ZO,ZO,ZO,ZO,ZO;
                                PB,NS,PS,PS,PS,PS,PB;
                                PB,PM,PM,PM,PS,PS,PB];
        end

        % 计算输入的隶属度及其在模糊子集中的位置
        function get_grad_membership(obj, erro, erro_c)
            if erro > obj.e_membership_values(1) && erro <= obj.e_membership_values(end)
                for i = 1:length(obj.e_membership_values)-1
                    if erro > obj.e_membership_values(i) && erro <= obj.e_membership_values(i+1)
                        obj.e_gradmembership(1) = -(erro-obj.e_membership_values(i+1))/(obj.e_membership_values(i+1)-obj.e_membership_values(i));
                        obj.e_gradmembership(2) = 1+(erro-obj.e_membership_values(i+1))/(obj.e_membership_values(i+1)-obj.e_membership_values(i));
                        obj.e_grad_idx(1) = i;
                        obj.e_grad_idx(2) = i + 1;
                    end
                end
            else
                if erro <= obj.e_membership_values(1)
                    obj.e_gradmembership(1) = 1;
                    obj.e_grad_idx(1) = 1;
                elseif erro > obj.e_membership_values(end)
                    obj.e_gradmembership(1) = 1;
                    obj.e_grad_idx(1) = 7;
                end
            end
            if erro_c > obj.ec_membership_values(1) && erro_c <= obj.ec_membership_values(end)
                for i = 1:length(obj.ec_membership_values)-1
                    if erro_c > obj.ec_membership_values(i) && erro_c <= obj.ec_membership_values(i+1)
                        obj.ec_gradmembership(1) = -(erro_c- obj.ec_membership_values(i+1)) / (obj.ec_membership_values(i+1) - obj.ec_membership_values(i));
                        obj.ec_gradmembership(2) = 1+(erro_c - obj.ec_membership_values(i+1))/(obj.ec_membership_values(i + 1) - obj.ec_membership_values(i));
                        obj.ec_grad_idx(1) = i;
                        obj.ec_grad_idx(2) = i + 1;
                    end
                end
            else
                if erro_c <= obj.ec_membership_values(1)
                    obj.ec_gradmembership(1) = 1;
                    obj.ec_grad_idx(1) = 1;
                elseif erro_c > obj.ec_membership_values(end)
                    obj.ec_gradmembership(1) = 1;
                    obj.ec_grad_idx(1) = 7;
                end
            end
        end

        % 计算pid参数增量总隶属度
        function getsumgrad(obj)
            for i=1:7
                obj.delta_kp_grad_sums(i) = 0;
                obj.delta_ki_grad_sums(i) = 0;
                obj.delta_kd_grad_sums(i) = 0;
            end
            for i=1:length(obj.e_grad_idx)
                for j=1:length(obj.ec_grad_idx)
                    idx_kp = obj.kp_rule_list(obj.e_grad_idx(i),obj.ec_grad_idx(j)) + 4;
                    idx_ki = obj.ki_rule_list(obj.e_grad_idx(i),obj.ec_grad_idx(j)) + 4;
                    idx_kd = obj.kd_rule_list(obj.e_grad_idx(i),obj.ec_grad_idx(j)) + 4;
                    obj.delta_kp_grad_sums(idx_kp) = obj.delta_kp_grad_sums(idx_kp) + obj.e_gradmembership(i)*obj.ec_gradmembership(j);
                    obj.delta_ki_grad_sums(idx_ki) = obj.delta_ki_grad_sums(idx_ki) + obj.e_gradmembership(i)*obj.ec_gradmembership(j);
                    obj.delta_kd_grad_sums(idx_kd) = obj.delta_kd_grad_sums(idx_kd) + obj.e_gradmembership(i)*obj.ec_gradmembership(j);
                end
            end
        end

        % 获取pid参数增量论域值
        function getout(obj)
            obj.qdelta_kp = 0;
            obj.qdelta_ki = 0;
            obj.qdelta_kd = 0;
            for i=1:7
                obj.qdelta_kp = obj.qdelta_kp + obj.kp_menbership_values(i)*obj.delta_kp_grad_sums(i);
                obj.qdelta_ki = obj.qdelta_ki + obj.ki_menbership_values(i)*obj.delta_ki_grad_sums(i);
                obj.qdelta_kd = obj.qdelta_kd + obj.kd_menbership_values(i)*obj.delta_kd_grad_sums(i);
            end
        end

        % 量化函数
        function y = quantization(~,maxvalue,minvalue,x)
            y = 6*(x-maxvalue)/(maxvalue-minvalue) + 3;
        end

        % 反量化函数
        function y = inverse_quantization(~,maxvalue,minvalue,qvalues)
            y = (maxvalue-minvalue)*(qvalues-3)/6 + maxvalue;
        end

        % 模糊pid控制主函数
        function [pidout,kpout,kiout,kdout,erro_c] = fuzzypidcontroller(obj,erro)
            erro_c = erro - obj.pre_err;
            obj.qerro = obj.quantization(obj.e_max,obj.e_min,erro);
            obj.qerro_c = obj.quantization(obj.ec_max,obj.ec_min,erro_c);
            obj.get_grad_membership(obj.qerro,obj.qerro_c)
            obj.getsumgrad()
            obj.getout()

            obj.delta_kp = obj.inverse_quantization(obj.delta_kpmax,obj.delta_kpmin,obj.qdelta_kp);
            obj.delta_ki = obj.inverse_quantization(obj.delta_kimax,obj.delta_kimin,obj.qdelta_ki);
            obj.delta_kd = obj.inverse_quantization(obj.delta_kdmax,obj.delta_kdmin,obj.qdelta_kd);

            obj.kp = obj.kp + obj.delta_kp;
            obj.ki = obj.ki + obj.delta_ki;
            obj.kd = obj.kd + obj.delta_kd;

            if obj.kp < 0
                obj.kp = 0;
            end
            if obj.ki < 0
                obj.ki = 0;
            end
            if obj.kd < 0
                obj.kd = 0;
            end
            kpout = obj.kp;
            kiout = obj.ki;
            kdout = obj.kd;

            pidout = obj.kp*erro + obj.ki*obj.erri + obj.kd*erro_c/obj.dt;
            obj.erri = obj.erri + erro*obj.dt;
            obj.pre_err = erro;

            if abs(pidout) > obj.threshold
                pidout = sign(pidout)*obj.threshold;
            end
        end
    end
end

main.m

% created by hyacinth on 2024/5/31
clc
clear
close all

%%
fs = 2e3;
dt = 1/fs;

pnums = 133;
pdens = [1,25,0];
P = tf(pnums,pdens);
dP = c2d(P,dt,"t");
[pnumz,pdenz] = tfdata(dP,"v");

tend = 5;
N = (tend/dt);
t = (1:N)*dt;
amp = 1000;
yd = amp*sin(2*pi*1*t);
% yd = amp*ones(1,N);

 kp0 = 0;
ki0 = 0;
kd0 = 0;

e_max = 0.05*amp;   % needs to be adjusted, maybe
e_min = -e_max;

ec_max = 1/fs*amp;  % needs to be adjusted, maybe
ec_min = -ec_max;

multi = 5e-4;       % needs to be adjusted, maybe
delta_kpmax = multi*1e3;
delta_kpmin = multi*-1e3;
delta_kimax = multi*1e0;
delta_kimin = multi*-1e0;
delta_kdmax = multi*1e-3;
delta_kdmin = multi*-1e-3;
threshold = 5000;

fuzzypid = FuzzyPID(e_max,e_min,ec_max,ec_min,kp0,ki0,kd0, ...
    delta_kpmax,delta_kpmin,delta_kimax,delta_kimin,delta_kdmax,delta_kdmin,threshold,dt);
zip = zeros(1,length(pdenz)-1);
u0 = 0;
for k = 1:N
    [y(k),zip] = filter(pnumz,pdenz,u0,zip);
    erro(k) = yd(k) - y(k);
    [u(k),kp(k),ki(k),kd(k),erro_c(k)] = fuzzypid.fuzzypidcontroller(erro(k));
    u0 = u(k);
end

figure
subplot(6,3,[1,4,7]);hold on;grid on
plot(t,yd,".-")
plot(t,y,".-")
legend("yd","y")
subplot(6,3,[10,13,16]);plot(t,u,".-");grid on;legend("u");
subplot(6,3,[2,5,8]);plot(t,erro,".-");grid on;legend("erro")
subplot(6,3,[11,14,17]);plot(t(2:end),erro_c(2:end),".-");grid on;legend("erro\_c")
subplot(6,3,[3,6]);plot(t,kp,".-");grid on;legend("kp")
subplot(6,3,[9,12]);plot(t,ki,".-");grid on;legend("ki")
subplot(6,3,[15,18]);plot(t,kd,".-");grid on;legend("kd")