Дифференциальная часть PID и его фильтр низких частот - область времени - PullRequest
Новая
       36

Дифференциальная часть PID и его фильтр низких частот - область времени

1 голос
/ 04 апреля 2019

Я запрограммировал PID в MATLAB: 

classdef PID < handle
    properties
        Kp = 0
        Ki = 0
        Kd = 0
        SetPoint = 1
        Dt = 0.01
    end

    properties (Access = private)
        IState = 0
        PreErr = 0
    end

    methods
        function obj = PID(Kp, Ki, Kd, SetPoint, Dt)
            if nargin == 0
                return;
            end
            obj.Kp = Kp;
            obj.Ki = Ki;
            obj.Kd = Kd;
            obj.SetPoint = SetPoint;
            obj.Dt = Dt;
        end

        function output = update(obj, measuredValue, t)
            err = obj.SetPoint - measuredValue;
            P = obj.getP(err);
            I = obj.getI(err);
            val = lowPass(obj,t);
            D = obj.getD(err*val);
            output = P + I + D;
        end

        function val = getP(obj, err)
            val = obj.Kp*err;
        end

        function val = getI(obj, err)
            obj.IState = obj.IState + err * obj.Dt;
            val = obj.Ki * obj.IState;
        end

        function val = getD(obj, err)
            val = obj.Kd * (err - obj.PreErr) / obj.Dt;            
            obj.PreErr = err;
        end

        function val = lowPass(obj,t)
            N = 10;
            val = 1-exp(-N*t);
        end
    end
end

И проверил его, используя случайный фильтр нижних частот в качестве установки: 

function r = getResponse(t)
r = 1 - exp(-5*t);
end

Код теста: 

sr = 1e2; % sampling rate 100Hz
st = 10; % sampling time 10s
ss = st*sr+1; % sample size
t = 0:1/sr:st; % time

input = ones(1,ss)*100;
output = zeros(1,ss);
measured = 0;

pid = PID(0,1,1,input(1),t(2)-t(1));
for i = 2:ss
    rPID(i) = pid.update(measured, t(i));
    output(i) = rPID(i)*getResponse(t(i));    
    measured = output(i);
end
figure
plot(t,output)
hold on;
plot(t,input)
plot(t,rPID)
legend('Output','Input','PID')

Обратите внимание, что для параметров установлено значение kp=0;ki=1;kd=1;.Я только тестирую дифференциальную часть здесь.Результат очень неправильный:

enter image description here

Обратите внимание, что ось Y масштабируется на 10 ^ 307.Он становится слишком большим, что через ~ 1.6 с значение ПИД превышает диапазон двойной точности и, следовательно, кривая останавливается.

Я убедился, что обе части P и I работают достаточно хорошо (см. этот вопрос, который я задал некоторое время назад ).

Из кривой для компонента D (см.рисунок ниже), хорошо видно, что он начинает сильно колебаться с самого начала;его значение достигает> 50k после 5-й отметки времени при 0,04 с:

enter image description here

Я почти уверен, что допустил ошибку при реализации нижнего проходафильтр, но я также заметил, что даже при удаленном фильтре нижних частот, дифференциальные значения по-прежнему ведут себя одинаково.

Чтобы иметь какую-то ссылку и сравнение, я также сделал симуляцию Simulink той же системыс использованием точно таких же коэффициентов усиления ПИД (т. е. kp=0;ki=1;kd=1;).Ниже приведена блок-схема (слева), рисунок для входа и выхода (вверху справа) и рисунок для значений ПИД (внизу справа)

enter image description here

Примечаниечто в блоках усиления нет верхнего / нижнего предела, а начальные входы / выходы установлены в нули.

Эти коэффициенты усиления ПИД далеко не оптимизированы, но они дают совершенно разные результаты в симуляции и кодированном ПИД.

Поэтому большой вопрос, я что-то здесь не так делаю?Почему разница между этими двумя результатами?

1 Ответ

2 голосов
/ 04 апреля 2019

Реализация фильтра нижних частот неверна. Разностное уравнение фильтра нижних частот выглядит следующим образом:

eq1: LowPass

Вызов функции getResponse может выглядеть так: 

pid = PID(0,1,1,input(1),t(2)-t(1)); 
for i = 2:ss
    rPID(i) = pid.update(measured, t(i));   
    alpha = getResponse(0.25,0.01);
    output(i) = rPID(i)*alpha+(1-alpha)*output(i-1);   
    measured = output(i);
end

Таким образом, getResponse эквивалентен альфа 

function r = getResponse(wc,Ts)
    r = 1 - exp(-wc*Ts);
end

Далее вам нужно изменить функцию lowPass в классе PID. 

    function output = update(obj, measuredValue)
        err = obj.SetPoint - measuredValue;
        P = obj.getP(err);
        I = obj.getI(err);
        val = lowPass(obj,err,0.1,0.01);
        D = obj.getD(val);
        output = P + I + D;
    end
    % ...
    function val = lowPass(obj,err,wc,Ts)
        alpha = getResponse(wc,Ts);
        output = err*alpha+(1-alpha)*obj.output_1;  
        obj.output_1 = output;
        val = output;
    end

The resulting plot:

...