Modelling of pipeline pressure control using a fuzzy controller

This example demonstrates the modelling of pipeline pressure control using a fuzzy controller.

The model file contains two hydraulic schemes with control systems. The upper circuit is controlled by a PID controller, its description can be found in the Community: Modelling Pipeline Pressure Control and Pipeline Pressure Control. The setpoint signal in both schemes is represented by a block Chart

The lower scheme has similar physical block parameters.

Schematic of the model:

The fuzzy controller is described using the block Engee Function. It allows you to use code in the Julia programming language. More details about working with Engee Function can be found in the corresponding documentation article.

To create a fuzzy controller, we will use the Julia programming language library called FuzzyLogic. To install it in Engee we need to execute the following code cell:

Pkg.add(["FuzzyLogic"])

Pkg.add("FuzzyLogic")

Run the installed FuzzyLogic library and the Plots plotting library.

using FuzzyLogic
using Plots

In the code cell below the construction of Mamdani fuzzy inference system is performed using the macro @mamfis. The construction of such a system includes the definition of membership functions. A function tipper is declared which takes one argument delta_p. This variable will represent the input data that comes from the adder that calculates the difference between the set pressure and the measured pressure.

The membership functions for pressure difference ranges and valve cross-sectional area outputs include:

• Gaussian membership function GaussianMF(x, y) - where x is the mean value defining the centre of the spread and y is the standard deviation
• triangular membership function TriangularMF(x, y, z), where x is the left corner of the triangle, z is the right corner and y is the peak of the triangle
• trapezoidal membership function TrapezoidalMF(x, y, z, k), where x is the lower left corner of the trapezoid, k is the lower right corner, y and z are the upper left and upper right corners respectively

The variable domain defines the ranges in which the membership functions will be defined.

Then fuzzy logic rules are defined, which link input and output variables. For example:

If delta_p is very high, then control_signal will be fully open. If delta_p is high, then control_signal will be partially open. And so on for the rest of the conditions.

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)

Plotting the membership functions for the input data:

plot(fis, :delta_p)

Plotting the membership functions for the output data of the data:

plot(fis, :control_signal)

In order to realise this fuzzy inference system as a control system, it is necessary to generate from it a function that is independent of the library and defined only by the basic constructs of the Julia language. For this purpose we can use the function compilefis(fis), where the argument fis is the fuzzy inference system defined earlier and written to a variable. The result of code generation with the help of this function is written to the variable 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)

The result written to the variable asd is overwritten in the code cell:

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)

Having defined the function tipper(), describing the fuzzy output system, we can plot the record of output data from input data, for this purpose we define the vector of pressure values from -50000 to 50000 Pa, as it was written earlier in the variable domain and pass it through the function:

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

Plotting the response curve of the resulting function:

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

We obtain a graph of the valve opening area dependence on the pressure difference between the set point and the sensor signal.

The function obtained with the help of code generation can be written to a file. To do this, the variable asd, which has the format Expr, is converted to String format for further writing:

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"

Writing a function in string format to a file:

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

We place the resulting file, of format .jl, with the function tipper(), into the block Engee Function and execute it using the function include(). To do this, go to the block parameters, click on the "Look under the mask" button, select the "Main" tab and click on "Edit source code", in the opened window we write:

struct Block <: AbstractCausalComponent; end

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

Next, the model is run and the data is plotted.

Defining the function to load and run the model:

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)

Running the simulation

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

Selecting data from the simout variable and writing them to variables:

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")

Writing the signals from the setpoint adjuster and from the pressure sensors of both circuits into variables:

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

Visualisation of simulation results

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

Having analysed the graph we can notice that unlike PID-regulator the fuzzy controller has much smaller overshoot and also shorter duration of transient process. But the fuzzy controller has a disadvantage, its signal can often be shifted from the set point signal by some value, but this can be corrected by more accurate adjustment of the belonging functions and solving rules.

Conclusion:

This example demonstrated the modelling of two physical objects with different variants of automatic control systems, in particular a PID controller and a fuzzy controller. A fuzzy inference system was constructed which, then, was converted into a function suitable for embedding in an Engee Function block. The result is a fuzzy controller that is as good as a PID controller.

Blocks used in example

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