Engee documentation
Notebook

Identification of a control system object with a fuzzy PID controller

In this example, two tasks are solved at once. First, we identify the control object based on experimental data. Secondly, we are implementing a fuzzy output system for the PID controller.

In [ ]:
Pkg.add(["ControlSystemIdentification", "FuzzyLogic"])
using ControlSystemIdentification, FuzzyLogic

Let me remind you that the identification workflow looks like this:

  1. Collecting experimental data

  2. Choosing the model structure

  3. Identification of the model

  4. Validation of the model

Identification of the management object

Collecting experimental data from the management facility

The simulation of the experimental installation of the control facility is presented in the plant.engee model. The control object is a discrete system with one input and one output. Noise and a slight delay have also been added to the model.

The identification of the object will be carried out by a transitional process, for this we will provide a step at the entrance of the system.

Plant--1753359832841.png

Let's record the identification data.:

In [ ]:
mode = 1; # Запись данных для идентификации
if "Plant" in [m.name for m in engee.get_all_models()]
    m = engee.open( "Plant" ) # загрузка модели
else
    m = engee.load( "Plant.engee" )
end
results = engee.run(m, verbose=true)
Building...
Progress 0%
Progress 100%
Progress 100%
Out[0]:
SimulationResult(
    "out" => WorkspaceArray{Float64}("Plant/out")
,
    "in" => WorkspaceArray{Float64}("Plant/in")

)
In [ ]:
dataout = collect(results["out"]);
datain = collect(results["in"]);
# Частота дискретизации системы из вектора времени
Ts = dataout.time[2] - dataout.time[1]
ident_data = iddata(dataout.value, datain.value, Ts)
plot(ident_data)
Out[0]:

Let's record the data for model validation.:

In [ ]:
mode = 0; # Запись данных для валидации
if "Plant" in [m.name for m in engee.get_all_models()]
    m = engee.open( "Plant" ) # загрузка модели
else
    m = engee.load( "Plant.engee" )
end
results = engee.run(m, verbose=true)
Building...
Progress 0%
Progress 100%
Progress 100%
Out[0]:
SimulationResult(
    "out" => WorkspaceArray{Float64}("Plant/out")
,
    "in" => WorkspaceArray{Float64}("Plant/in")

)
In [ ]:
dataout = collect(results["out"]);
datain = collect(results["in"]);
# Частота дискретизации системы из вектора времени
Ts = dataout.time[2] - dataout.time[1]
val_data = iddata(dataout.value, datain.value, Ts)
plot(val_data)
Out[0]:

Model structure selection and identification

We will identify a linear stationary system in the state space. Let's choose the order of the system equal to 3.

For identification, we will choose the PEM (prediction-error method) prediction error method. This is a simple but effective algorithm for identifying linear stationary systems with discrete time in the form of a state space. The method also allows you to get the initial states of the system.

In [ ]:
nx = 3; # порядок системы
sys, x0 = newpem(ident_data,nx)
Iter     Function value   Gradient norm 
     0     3.784040e-04     9.789230e-03
 * time: 0.19949603080749512
    50     2.082491e-04     4.272067e-04
 * time: 0.2358410358428955
   100     2.077983e-04     4.445332e-03
 * time: 0.3109588623046875
   150     2.068965e-04     6.060038e-03
 * time: 0.39757585525512695
   200     2.058789e-04     3.800742e-03
 * time: 0.4177060127258301
   250     2.057405e-04     2.125884e-03
 * time: 0.45604991912841797
   300     2.056074e-04     6.434034e-04
 * time: 0.48035192489624023
   350     2.055462e-04     3.334684e-04
 * time: 0.5061249732971191
   400     2.053596e-04     1.005185e-03
 * time: 0.5335268974304199
   450     2.051661e-04     9.945758e-03
 * time: 0.5657720565795898
   500     2.050343e-04     1.161464e-03
 * time: 0.6026959419250488
   550     2.048684e-04     1.551165e-02
 * time: 0.6348810195922852
   600     2.047113e-04     1.148014e-02
 * time: 0.6805129051208496
   650     2.046216e-04     2.184856e-02
 * time: 0.6981239318847656
   700     2.045511e-04     8.645446e-03
 * time: 0.7165648937225342
   750     2.044672e-04     1.919704e-02
 * time: 0.7359099388122559
   800     2.043937e-04     2.137336e-02
 * time: 0.7552690505981445
   850     2.042809e-04     1.807003e-02
 * time: 0.7735428810119629
   900     2.042222e-04     2.871188e-02
 * time: 0.7964169979095459
   950     2.041250e-04     1.880932e-02
 * time: 0.8146660327911377
  1000     2.040802e-04     1.426929e-01
 * time: 0.8476660251617432
  1050     2.039567e-04     7.880180e-02
 * time: 0.8641018867492676
  1100     2.039079e-04     1.105847e-01
 * time: 0.8808779716491699
  1150     2.038122e-04     3.309354e-01
 * time: 0.89768385887146
  1200     2.037250e-04     5.343228e-01
 * time: 0.9144978523254395
  1250     2.036590e-04     2.938638e-01
 * time: 0.9343979358673096
  1300     2.036076e-04     1.312110e+00
 * time: 0.9605119228363037
  1350     2.035490e-04     1.029981e+00
 * time: 0.9860918521881104
  1400     2.034784e-04     1.900128e+00
 * time: 1.0114190578460693
  1450     2.034470e-04     1.202898e+00
 * time: 1.0391180515289307
  1500     2.034154e-04     1.747894e-01
 * time: 1.0554909706115723
  1550     2.033921e-04     1.671672e+00
 * time: 1.0719490051269531
  1600     2.033679e-04     3.414074e+00
 * time: 1.0886778831481934
  1650     2.033440e-04     9.629020e+00
 * time: 1.1050329208374023
  1700     2.033089e-04     6.344263e+00
 * time: 1.1232068538665771
  1750     2.032809e-04     3.688618e+00
 * time: 1.1415619850158691
  1800     2.032635e-04     1.999306e+01
 * time: 1.1600329875946045
Out[0]:
(sys = PredictionStateSpace{ControlSystemsBase.Discrete{Float64}, ControlSystemsBase.StateSpace{ControlSystemsBase.Discrete{Float64}, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}
A = 
 0.3140163495012443     3.7117022260861     0.0
 0.11023755964659158    1.6140225481406154  0.07598141692800656
 0.0                  -12.659066417167994   0.21882581128873133
B = 
 -3.320403895114945
  1.7883472038146635
 -7.720330696701567
C = 
 -6.177421822749285  -26.798969096921468  -3.5509173987959497
D = 
 0.0

Sample Time: 0.1 (seconds)
Discrete-time state-space model, x0 = [1.2762453764639579, -0.09723096657064308, -1.486442954748024], res =  * Status: success

 * Candidate solution
    Final objective value:     2.032499e-04

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 0.00e+00 ≤ 0.0e+00
    |x - x'|/|x'|          = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 1.0e-16
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 1.83e-02 ≰ 1.0e-12

 * Work counters
    Seconds run:   1  (vs limit 100)
    Iterations:    1841
    f(x) calls:    2387
    ∇f(x) calls:   1841
, nlp = Float64[])

Validation of the model

Let's check the identification on the same data that was used for identification.

In [ ]:
simplot(ident_data, sys)
Out[0]:

And we'll check it on a special validation dataset.

In [ ]:
simplot(val_data, sys)
Out[0]:

The method showed a good identification result.

Fuzzy PID controller control system

Let's return to modeling a control system with a fuzzy controller. We have obtained a model of the control object from experimental data.

Fuzzycontroller--1753369066574.png

The identification result will be used in the Discrete State Space block to model the control object.

Fuzzycontroller--1753369100219.png

It remains to add the controller to the control system circuit. This example uses the following fuzzy logic controller structure:

Fuzzycontroller--1753362744936.png

Normalized error values are received at the input of the fuzzy output system and the derived error . The input values are normalized using scaling factors. and so that they are in the range [-1,1].

The fuzzy output system outputs a result in the range [-1, 1], which is scaled by coefficients and .

In [ ]:
CE = 10;
CD = 4.9829;
C0 = 2.5179;
C1 = 1.8103;

The fuzzy inference system is represented as a block of a Lookup table. The process of forming the response surface for this table is shown below:

In [ ]:
fis = @sugfis function FCtrl(E, delE)::U
    E := begin
        domain = -10.0:10.0
        N = GaussianMF(-10.0, 7.0)
        P = GaussianMF(10.0, 7.0)
    end

    delE := begin
        domain = -10.0:10.0
        N = GaussianMF(-10.0, 7.0)
        P = GaussianMF(10.0, 7.0)
    end    

    U := begin
        domain = -20.0:20.0
        Min = -20.0
        Zero = 0.0
        Max  = 20.0
    end

   # Работающие правила
   E == N  && delE == N  -->  U == Min
   E == N  && delE == P  -->  U == Zero
   E == P  && delE == N  -->  U == Zero
   E == P  && delE == P  -->  U == Max

end;
In [ ]:
plot(fis, :E, ylabel="Функция принадлежности", w = 2)
Out[0]:
In [ ]:
plot(fis, :delE, ylabel="Функция принадлежности", w = 2)
Out[0]:

Let's create a response surface in order to place it in the Lookup table block.

In [ ]:
errBrkpts = -10.0:0.5:10.0 
rateBrkpts = -10.0:0.5:10.0
LookUpTableData = zeros(length(errBrkpts), length(rateBrkpts))
i = 0;
for e in errBrkpts
    i = i+1
    j = 0;
    for de in rateBrkpts
        j = j+1
        LookUpTableData[i, j] = fis([e, de])[:U]
    end
end
surface(errBrkpts, rateBrkpts, LookUpTableData, xlabel="ΔE", ylabel="E", zlabel="U", title="Поверхность отклика", camera=(30, 30))
Out[0]:

Let's run the control system model and check the operation of the regulator.

In [ ]:
if "Fuzzycontroller" in [m.name for m in engee.get_all_models()]
    m = engee.open( "Fuzzycontroller" ) # загрузка модели
else
    m = engee.load( "Fuzzycontroller.engee" )
end
fuzzyresults = engee.run(m, verbose=true)
Building...
Progress 0%
Progress 71%
Progress 100%
Progress 100%
Out[0]:
SimulationResult(
    "out" => WorkspaceArray{Float64}("Fuzzycontroller/out")
,
    "in" => WorkspaceArray{Float64}("Fuzzycontroller/in")

)
In [ ]:
fuzzydataout = collect(fuzzyresults["out"]);
fuzzydatain = collect(fuzzyresults["in"]);

In the graphs, we see the response of the control system to the stepwise impact.

In [ ]:
plot(fuzzydataout.time, fuzzydataout.value, titel = "Система управления с нечетким ПИД-регулятором", label="выход")
plot!(fuzzydatain.time, fuzzydatain.value, label="вход")
Out[0]:

Conclusion

As a result, we got a working control system with a fuzzy PID controller. The model of the control object was not known in advance, but was obtained by identification based on experimental data.