Моделирование управления давлением в трубопроводе с помощью нечеткого регулятора

Открыть пример в Engee

В данном примере продемонстрировано моделирование управления давлением в трубопроводе с помощью нечёткого регулятора.

В файле модели содержатся две гидравлические схемы с системами управления. Верхняя схема управляется PID-регулятором, её описание можно найти в Сообществе: Моделирование управления давлением в трубопроводе и Регулирование давления в трубопроводе. Сигнал задатчика в обеих схемах представлен блоком Chart

Нижняя схема имеет аналогичные параметры физических блоков.

Схема модели:

Нечёткий регулятор описан с помощью блока Engee Function. Он позволяет использовать код на языке программирования Julia. Подробнее о работе с Engee Function можно узнать в соответствующей статье документации.

Чтобы создать нечёткий регулятор воспользуемся библиотекой языка программирования Julia, которая называется FuzzyLogic. Для её установки в Engee требуется выполнить следующую кодовую ячейку:

Pkg.add(["FuzzyLogic"])

Pkg.add("FuzzyLogic")

Запускаем установленную библиотеку FuzzyLogic и библиотеку для построения графиков Plots.

using FuzzyLogic
using Plots

В кодовой ячейке ниже осуществляется построение системы нечёткого вывода Мамдани (Mamdani fuzzy inference system) с помощью макроса @mamfis. Построение такой системы включает в себя определение функций принадлежности (Membership functions). Объявляется функция tipper, которая принимает один аргумент delta_p. Эта переменная будет представляет собой входные данные, которые поступают с сумматора, вычисляющего разность между заданным и измеренным давлением.

Функции принадлежности для диапазонов разностей давлений и выходных сигналов площади сечения клапана включают в себя:

• гауссовскую функцию принадлежности GaussianMF(x, y) - где x - среднее значение, определяющее центр распраделения, а y - стандартное отклонение
• треугольную функцию принадлежности TriangularMF(x, y, z), где x - левый угол треугольника, z - правый, а y - пик треугольника
• трапециевидную функцию принадлежности TrapezoidalMF(x, y, z, k), где x - левый нижний угол трапеции, k - правый нижний угол, y и z - левый верхний и правый верхний углы соответственно

В переменной domain определяются диапазоны в которых будут определены функции принадлежности.

Далее определяются правила нечёткой логики, которые связывают входные и выходные переменные. Например:

Если delta_p очень высокое, то control_signal будет полностью открытым. Если delta_p высокое, то control_signal будет частично открытым. И так далее для остальных условий.

fis = @mamfis function tipper(delta_p)::control_signal
    delta_p := begin
        domain = -50000.0:50000.0
        very_high = GaussianMF(17000.0, 3000.0)  # Очень высокое давление
        high = GaussianMF(15000.0, 4120.0)       # Высокое давление
        moderate = GaussianMF(-6000.0, 3000.0)       # Нормальное давление
        low = GaussianMF(-15000.0, 3000.0)       # Низкое давление
        very_low = GaussianMF(-30000.0, 3000.0)  # Очень низкое давление
    end

    control_signal := begin
        domain = 0.0:0.005
        fully_closed = TriangularMF(0.00000000001, 0.00000001, 0.000001)  # Полностью закрытая заслонка
        partially_closed = TriangularMF(0.0000009, 0.000005, 0.00001)  # Частично закрытая заслонка
        neutral = TriangularMF(0.00001, 0.00002, 0.0003)  # Нейтральное положение (половина открыта)
        partially_open = TrapezoidalMF(0.0003, 0.0009, 0.003, 0.0044)  # Частично открытая заслонка
        fully_open = TriangularMF(0.004, 0.0047, 0.005)  # Полностью открытая заслонка
    end

    # Работающие правила
    delta_p == very_high --> control_signal == fully_open
    delta_p == high --> control_signal == partially_open
    delta_p == moderate --> control_signal == neutral
    delta_p == low --> control_signal == partially_closed
    delta_p == very_low --> control_signal == fully_closed
end

tipper

Inputs:
-------
delta_p ∈ [-50000.0, 50000.0] with membership functions:
    very_high = GaussianMF{Float64}(17000.0, 3000.0)
    high = GaussianMF{Float64}(15000.0, 4120.0)
    moderate = GaussianMF{Float64}(-6000.0, 3000.0)
    low = GaussianMF{Float64}(-15000.0, 3000.0)
    very_low = GaussianMF{Float64}(-30000.0, 3000.0)


Outputs:
--------
control_signal ∈ [0.0, 0.005] with membership functions:
    fully_closed = TriangularMF{Float64}(1.0e-11, 1.0e-8, 1.0e-6)
    partially_closed = TriangularMF{Float64}(9.0e-7, 5.0e-6, 1.0e-5)
    neutral = TriangularMF{Float64}(1.0e-5, 2.0e-5, 0.0003)
    partially_open = TrapezoidalMF{Float64}(0.0003, 0.0009, 0.003, 0.0044)
    fully_open = TriangularMF{Float64}(0.004, 0.0047, 0.005)


Inference rules:
----------------
delta_p is very_high --> control_signal is fully_open
delta_p is high --> control_signal is partially_open
delta_p is moderate --> control_signal is neutral
delta_p is low --> control_signal is partially_closed
delta_p is very_low --> control_signal is fully_closed


Settings:
---------
- MinAnd()
- MaxOr()
- MinImplication()
- MaxAggregator()
- CentroidDefuzzifier(100)

Вывод на график функций принадлежности для входных данных:

plot(fis, :delta_p)

Вывод на график функций принадлежности для выходных данных данных:

plot(fis, :control_signal)

Для того чтобы реализовать данную систему нечёткого вывода в виде системы управления, необходимо сгенерировать из неё функцию, не зависящую от библиотеки и определяемую только базовыми конструкциями языка Julia. Для этого можно использовать функцию compilefis(fis), где аргумент fis это система нечёткого вывода, определённая ранее и записанная в переменную. Результат кодогенерации с помощью данной функции записывается в переменную asd:

asd = compilefis(fis)

:(function tipper(delta_p)
      very_high = exp(-((delta_p - 17000.0) ^ 2) / 1.8e7)
      high = exp(-((delta_p - 15000.0) ^ 2) / 3.39488e7)
      moderate = exp(-((delta_p - -6000.0) ^ 2) / 1.8e7)
      low = exp(-((delta_p - -15000.0) ^ 2) / 1.8e7)
      very_low = exp(-((delta_p - -30000.0) ^ 2) / 1.8e7)
      ant1 = very_high
      ant2 = high
      ant3 = moderate
      ant4 = low
      ant5 = very_low
      control_signal_agg = collect(LinRange{Float64}(0.0, 0.005, 101))
      @inbounds for (i, x) = enumerate(control_signal_agg)
              fully_closed = max(min((x - 1.0e-11) / 9.99e-9, (1.0e-6 - x) / 9.9e-7), 0)
              partially_closed = max(min((x - 9.0e-7) / 4.1000000000000006e-6, (1.0e-5 - x) / 5.0e-6), 0)
              neutral = max(min((x - 1.0e-5) / 1.0e-5, (0.0003 - x) / 0.00028), 0)
              partially_open = max(min((x - 0.0003) / 0.0006000000000000001, 1, (0.0044 - x) / 0.0014000000000000002), 0)
              fully_open = max(min((x - 0.004) / 0.0007000000000000001, (0.005 - x) / 0.0002999999999999999), 0)
              control_signal_agg[i] = max(max(max(max(min(ant1, fully_open), min(ant2, partially_open)), min(ant3, neutral)), min(ant4, partially_closed)), min(ant5, fully_closed))
          end
      control_signal = ((2 * sum((mfi * xi for (mfi, xi) = zip(control_signal_agg, LinRange{Float64}(0.0, 0.005, 101)))) - first(control_signal_agg) * 0.0) - last(control_signal_agg) * 0.005) / ((2 * sum(control_signal_agg) - first(control_signal_agg)) - last(control_signal_agg))
      return control_signal
  end)

Результат, записанный в переменную asd переписываем в кодовую ячейку:

function tipper(delta_p)
    very_high = exp(-((delta_p - 17000.0) ^ 2) / 1.8e7)
    high = exp(-((delta_p - 15000.0) ^ 2) / 3.39488e7)
    moderate = exp(-((delta_p - -6000.0) ^ 2) / 1.8e7)
    low = exp(-((delta_p - -15000.0) ^ 2) / 1.8e7)
    very_low = exp(-((delta_p - -30000.0) ^ 2) / 1.8e7)
    ant1 = very_high
    ant2 = high
    ant3 = moderate
    ant4 = low
    ant5 = very_low
    control_signal_agg = collect(LinRange{Float64}(0.0, 0.005, 101))
    @inbounds for (i, x) = enumerate(control_signal_agg)
            fully_closed = max(min((x - 1.0e-11) / 9.99e-9, (1.0e-6 - x) / 9.9e-7), 0)
            partially_closed = max(min((x - 9.0e-7) / 4.1000000000000006e-6, (1.0e-5 - x) / 5.0e-6), 0)
            neutral = max(min((x - 1.0e-5) / 1.0e-5, (0.0003 - x) / 0.00028), 0)
            partially_open = max(min((x - 0.0003) / 0.0006000000000000001, 1, (0.0044 - x) / 0.0014000000000000002), 0)
            fully_open = max(min((x - 0.004) / 0.0007000000000000001, (0.005 - x) / 0.0002999999999999999), 0)
            control_signal_agg[i] = max(max(max(max(min(ant1, fully_open), min(ant2, partially_open)), min(ant3, neutral)), min(ant4, partially_closed)), min(ant5, fully_closed))
        end
    control_signal = ((2 * sum((mfi * xi for (mfi, xi) = zip(control_signal_agg, LinRange{Float64}(0.0, 0.005, 101)))) - first(control_signal_agg) * 0.0) - last(control_signal_agg) * 0.005) / ((2 * sum(control_signal_agg) - first(control_signal_agg)) - last(control_signal_agg))
    return control_signal
end

tipper (generic function with 2 methods)

Определив функцию tipper(), описывающую систему нечёткого вывода, можем построить график записимости выходных данных от входных, для этого определяем вектор значений давления от -50000 до 50000 Па, как было записано ранее в переменной domain и пропускаем его через функцию:

x = collect(range(-50000.0,50000.0,10001));

Построение кривой отклика получившейся функциии:

plot(x, tipper.(x), linewidth=3)

Получем график зависимости площади открытия клапана от разности давлений между задатчиком и сигналом с датчика.

Функцию, полученную с помощью кодогенерации, можно записать в файл. Для этого переменная asd имеющая формат Expr, конвертируется в формат String, для дальнейшей записи:

text_function = string(asd)

"function tipper(delta_p)\n    very_high = exp(-((delta_p - 17000.0) ^ 2) / 1.8e7)\n    high = exp(-((delta_p - 15000.0) ^ 2) / 3.39488e7)\n    moderate = exp(-((delta_p - -6000.0) ^ 2) / 1.8e7)\n    low = exp(-((delta_p - -15000.0) ^ 2) / 1.8e7)\n    very_low = exp(-((delta_p" ⋯ 977 bytes ⋯ "xi for (mfi, xi) = zip(control_signal_agg, LinRange{Float64}(0.0, 0.005, 101)))) - first(control_signal_agg) * 0.0) - last(control_signal_agg) * 0.005) / ((2 * sum(control_signal_agg) - first(control_signal_agg)) - last(control_signal_agg))\n    return control_signal\nend"

Запись функции в строковом формате в файл:

f = open("/user/start/examples/controls/liquid_pressure_regulator_fuzzy/liquid_pressure_regulator_fuzzy.jl","w") # создание файла
write(f, text_function) # запись в файл
close(f) # закрытие файла

Получившийся файл, формата .jl, с функцией tipper(), помещаем в блок Engee Function и выполняем с помощью функции include(). Чтобы это осуществить, необходимо зайти в параметры блока, нажать на кнопку "Посмотреть под маску", выбрать вкладку "Main" и нажать "Редактировать исходный код", в открывшемся окне вписываем:

struct Block <: AbstractCausalComponent; end

function (c::Block)(t::Real, delta_p)
    include("liquid_pressure_regulator_fuzzy.jl")
    return control_signal
end

Далее осуществляется запуск модели и вывод данных на графики.

Определение функции для загрузки и запуска модели:

function start_model_engee()
    try
        engee.close("liquid_pressure_regulator_fuzzy", force=true) # закрытие модели 
        catch err # в случае, если нет модели, которую нужно закрыть и engee.close() не выполняется, то будет выполнена её загрузка после catch
            m = engee.load("$(@__DIR__)/liquid_pressure_regulator_fuzzy.engee") # загрузка модели
        end;

    try
        engee.run(m) # запуск модели
        catch err # в случае, если модель не загружена и engee.run() не выполняется, то будут выполнены две нижние строки после catch
            m = engee.load("$(@__DIR__)/liquid_pressure_regulator_fuzzy.engee") # загрузка модели
            engee.run(m) # запуск модели
        end
end

start_model_engee (generic function with 1 method)

Запуск симуляции

try
    start_model_engee() # запуск симуляции с помощью специальной функции, реализованной выше
    catch err
    end;

Выделение из переменной simout данных и их запись в переменные:

sleep(5)
result = simout;
res = collect(result)

4-element Vector{WorkspaceArray}:
 WorkspaceArray("liquid_pressure_regulator_fuzzy/Датчик давления.1")
 WorkspaceArray("liquid_pressure_regulator_fuzzy/Нечёткий регулятор.1")
 WorkspaceArray("liquid_pressure_regulator_fuzzy/Датчик давления-1.1")
 WorkspaceArray("liquid_pressure_regulator_fuzzy/Chart.RefPress")

Запись в переменные сигналов с задатчика и с датчиков давления обеих схем:

PID_регулятор = collect(res[1])
Нечёткий_регулятор = collect(res[3])
Сигнал_задатчика = collect(res[4]);

Визуализация результатов моделирования

plot(PID_регулятор[:,1], PID_регулятор[:,2], linewidth=3, label="PID-регулятор")
plot!(Нечёткий_регулятор[:,1], Нечёткий_регулятор[:,2], linewidth=3, label="Нечёткий регулятор")
plot!(Сигнал_задатчика[:,1], Сигнал_задатчика[:,2], linewidth=3, label="Сигнал задатчика")

Проанализировав график можно заметить, что в отличие от PID-регулятора у нечёткого регулятора значительно меньшее перерегулирование, а также меньше длительность переходного процесса. Но у нечёткого регулятора есть и минус, его сигнал, зачастую, может быть смещён от сигнала задатчика на некоторую величину, но это можно исправить более точной настройкой функций принадлежности и решающих правил.

Вывод:

В данном примере было продемонстрировано моделирование двух физических объектов с разными вариантами систем автоматического управления, в частности с PID-регулятором и с нечётким регулятором. Была построена система нечёткого вывода, которая, затем, была преобразована в функцию, подходящую для встраивания в блок Engee Function. В результате получился нечёткий регулятор, который по своим характеристикам не уступает PID-регулятору.

Блоки, использованные в примере

• Pressure Sensor (IL)
• PID Controller
• Delay
• Chart
• Engee Function
• Pipe (IL)